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1.  Introduction 


The  purpose  of  this  Grant  is  to  develop  theoretical  tools  for  performance 
assessment  of  Plasma  Propulsion  devices.  These  models  can  be  used  for  design  purposes, 
as  correlators  of  experimental  data,  and  as  pointers  to  necessary  technology  improvements. 

In  our  first  Yearly  Progress  Report  1^1  we  showed  results  of  a  simplified  arcjet 
model  which  was  then  being  finalized  by  Graduate  Student  Atsuko  Sakamoto.  The  final 
results  of  this  work  have  since  appeared  in  Ms.  Sakamoto’s  M.S.  Thesis  and  have  been 
presented  in  a  paper  at  the  29th  Joint  Propulsion  Conferencef^l.  The  method  is  capable  of 
good  (±20%)  a  priori  predictions  of  thrust  and  efficiency  over  a  wide  range  of  parameters, 
but  requires  some  judgment  in  the  selection  of  input  parameters,  and  provides  little  detailed 
information  on  the  internal  flow  parameters.  Also  in  the  first  Yearly  Report,  we  showed 
the  development  and  preliminary  results  of  a  detailed  arcjet  numerical  model  which  was  part 
of  Scott  Miller's  Ph.D.  Thesis  work.  This  model  is  now  essentially  complete,  and  results 
have  also  been  presented  at  the  Joint  Propulsion  Conferences'll.  This  is  probably  the  rr»st 
accurate  description  available  of  the  hydrogen  arcjet  physics,  and  provides  first-principles 
performance  predictions  as  well  as  internal  parameter  distributions  which  illuminate  the 
device's  operation  and  suggest  routes  for  improvement.  The  JPC  paper  summarizes  these 
results,  and  is  appended  for  further  detail.  Additional  results,  dealing  principally  with  the 
thermal  coupling  between  the  gas  and  the  arcjet  body,  are  still  being  generated,  and  will  be 
the  subject  of  a  subsequent  paper  to  be  presented  at  the  23rd  International  Electric 
Propulsion  Conference.  Scott  Miller's  Ph.D.  Thesis  is  expected  to  be  available  by  SepL. 
1993. 


In  addition  to  the  arcjet  work,  we  have  also  developed  during  this  2nd  Grant  Year  a 
one-dimensional  model  for  a  Hall  thruster,  which  has  yielded  excellent  results  when 
compared  to  experimental  data.  This  has  been  the  subject  of  a  M.S.  Thesis  by  Chris 
Lentzl^l,  and  also  of  a  paper  presented  at  the  Joint  Propulsion  Conferencef^J.  Although 
partial  attempts  had  been  made  before  at  modeling  Hall  thruster  operation,  it  appears  that 
ours  is  the  only  truly  quantitative  and  accurate  theory  in  existence  at  this  time.  The  work  is 
described  in  the  JPC  paper,  which  is  appended  to  this  Report,  and  in  more  detail  in  Ch. 
Lentz'  Thesis,  which  is  expected  in  August  1993. 

The  good  results  obtained  in  both  the  arcjet  and  the  Hall  thruster  tasks  have  opened 
the  way  for  further  refinements  and  for  explorations  of  performance-raising  modifications. 
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Thus,  work  is  now  being  initiated  (under  a  new  AFOSR  Grant)  on  mixed  N2-H2  arcjet 
modeling,  examination  of  alkali  seeding  in  arcjets,  and  2-D  Hall  Thruster  models. 

Although  the  Grant  being  reported  on  did  not  cover  work  on  MPD  thrusters,  work 
that  had  been  initiated  under  earlier  AFOSR  Grants  has  been  continued.  In  particular,  the 
2-D  numerical  model  of  MPD  which  was  developed  by  Eli  Niewood  to  examine  in  detail 
the  genesis  of  near-anode  voltage  drops  has  been  completed.  This  was  the  subject  of  Dr. 
Niewood's  Thesisl^I,  and  was  also  presented  at  the  Joint  Propulsion  Conferencel^J.  This 
work  constitutes  the  first  coherent  explanation  for  the  appearance  and  magnitude  of  the 
anode  drop  in  MPD  thrusters,  which  is  by  far  their  largest  power  loss  and  is  also  the  source 
of  difficult  anode  thermal  management  problems.  The  JPC  paper  is  appended  to  this 
Report,  and  a  copy  of  Dr.  Niewood's  Thesis  is  also  included. 
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2.  Numerical  Model  for  Hydrogen.  Arc iets 

Previous  numerical  modeling  of  arcjei  thrusters  has  focused  on  the  development  of 
1-D,  2-D,  and  axisymmetric  models  with  relatively  simple  physics  and  geometries.  The 
level  of  detail  has  ranged  from  1-D  models  to  coupled  quasi-analytic  models  of  the  inner 
(arc)  and  outer  (cold  gas)  flows  P.lO.l  l)^  to  simplified  axisymmetric  space-marching 
techniques  and  finally  to  2-D  and  axisymmetric  viscous  codes  which  begin  to 
incorporate  most  of  the  detailed  physical  processes  (13.14,15]  jhe  latest  research  has 
obtained  results  which  variously  include  ohmic  heating,  electron  heat  conduction,  thermal 
and  ionizational  nonequilibrium,  and  empirical  models  of  radiation  losses.  Issues  left  out 
of  these  models  include  viscous  and  diffusive  effects  (of  importance  in  these  very  low 
Reynolds  number  devices),  coupling  of  the  block  thermal  response  (which  affects  the  gas 
temperature  around  the  arc,  and  hence  the  specific  impulse),  and,  more  importantly, 
electron  thermal  non-equilibrium  {T^  *^t)  detailed  arc  attachment  physics.  If 

thermal  equilibrium  is  forced,  as  has  been  done  so  far,  the  electrical  conductivity  outside 
the  arc,  and  panicularly  near  the  anode  wall,  is  nearly  zero,  because  the  single  temperature 
is  held  down  by  cooling  to  the  wall.  This  makes  it  impossible  to  obtain  steady  solutions 
with  a  prescribed  current  of  the  right  magnitude,  unless  an  artificial  "floor"  is  introduced 
for  the  conductivity.  Unfortunately  neither  the  magnitude  nor  the  spatial  distribution  of  this 
artificial  can  be  gleaned  from  such  a  model,  and  this  throws  in  doubt  the  voltage 

prediction  as  well  as  that  of  the  anode  heat  flux  and  its  distribution.  To  some  extent,  this 
has  been  mitigated  by  anificially  imposing  (as  well)  the  current  distribution  on  the  anode, 
based  on  empirical  information. 

In  our  work,  we  have  decoupled  Te  from  Tg  by  separately  balancing  electron  and 
gas  energies.  The  results  show  that  Te  remains  close  to  Tg  in  the  arc,  where  the  collisional 
coupling  is  strong,  but  whereas  the  near- anode  Tg  is  of  the  order  of  1000-3000K,  the 
corresponding  Te  is  of  the  order  of  20,000K  there,  and,  in  fact,  that  a  continuous  "bridge" 
zone  of  high  Te  connects  the  arc  to  the  anode.  This  is  shown  clearly  in  Fig.  1,  which  is  a 
radial  cut  in  the  attachment  region.  Since  the  electrical  conductivity  is 
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Fig.  1;  Radial  Profiles  of  Electron  and  Heavy  Species  Temperatures  0.25nun  Downstream 

of  the  Constrictor  Exit 

governed  by  Te,  this  bridge  all  the  way  to  the  wall  allows  the  current  connection  to  be 
made,  with  a  net  voltage  drop  which  is  very  close  to  that  measured  experimental  (see  the 
attached  paper,  AIAA  93-2101,  for  details). 

Not  as  essential,  but  also  imponant  for  the  physics  of  the  attached  region  is  the 
allowance  for  ambipolar  diffusion  of  electron-ion  pairs  from  a  region  of  intense  net 
ionization  about  0.1mm  from  the  wall  to  the  wall  itself,  bridging  any  possible  electron-poor 
thin  layer  there.  This  is  shown  in  Fig.  2,  where  the  doirunant  terms  in  the  electron  density 
balance  are  profiled. 
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Figure  2-  Radial  Profiles  of  Some  Terms  in  the  Electron  Density  Equation  in  the  Current 
Attachment  Region  Just  Downstream  of  the  Constrictor  Exit 


When  these  effects  are  tracked  through  the  2-D  computational  domain,  the  it- 
attachment  distribution  in  the  region  downstream  of  the  constrictor  is  obtained  naturally  in  a 
manner  which  is  realistic,  although  no  detailed  experimental  distributions  are  available  for 
hot  anodes  (Sec  Fig.  4  of  the  attached  paper  AIAA  93-2101).  The  veracity  of  the  Te 
calculation  is  also  proven  by  comparing  the  predicted  Te  profile  a  the  thruster  exit  (Fig.  3) 
to  data  of  Haskins  et  al  and  of  Spores  In  both  cases  the  measured  Te  and  its  radial 
profile  agree  within  0. 1  eU  with  our  results;  the  data  of  Ref.  16  were  taken  a  few  mm. 
downstream  of  the  exit  plane,  which  may  account  for  most  of  the  difference. 

Despite  these  important  improvements  to  the  anode  attachment  problem,  no 
complete  success  can  yet  be  claimed.  On  one  hand,  we  do  slightly  under-predict  voltage 
(by  7V  out  of  1 12  in  the  nominal  case,  but  exploratory  calculations  for  other  conditions 
give  voltage  errors  of  up  to  19V.  in  some  instances).  More  importantly,  we  are  still  forced 
to  restrict  attachment  to  the  region  beyond  the  constrictor,  and  if  this  limitation  is  removed, 
attachment  gradually  creeps  upstream  without  apparent  limits.  The  mechanism  which 
anchors  the  arc  foot  is  not  completely  captured,  and  may  require  finer  detail  of  the  sheath 
and  the  surface  chemistry  to  be  modeled. 
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A  broad  measure  of  success  for  the  model  is  its  performance  predictive  ability.  As 
Fig.  1 1  of  the  included  Paper  AIAA  93-2101  shows,  the  thrust  (or  the  specific  impulse)  is 
within  5-10%  of  the  experimental  data  in  all  cases  tested,  consistently  showing  some 
overprediction. 

The  situation  for  voltage  is  similar  (Table  6  of  the  Paper).  A  wider  set  of  conditions 
is  currently  being  explored,  and  the  results  will  be  documented  in  an  EPC  paper  in 
September. 

The  other  improvement  contributed  by  our  model  is  the  ability  to  accurately  model 
the  thermal  response  of  a  radiatively  cooled  thruster  body,  and  to  couple  this  to  the  fluid 
model.  The  thermal  time  constant  is  of  the  order  of  seconds,  compared  to  the  -50  psec 
time  scale  for  flow  development,  and  this  makes  it  impossible  to  numerically  track  both 
time  evolutions  in  an  interactive  manner.  Instead,  the  calculation  is  iterative,  alternating 
between  a  flow  solution  with  a  frozen  wall  temperature  and  a  steady  state  thermal 
computation  using  as  input  the  heat  flux  from  the  latest  flow  update.  This  is  mainly 
imponant  in  setting  the  outer  gas  temperature,  which  controls  the  thruster  specific  impulse. 
In  fact,  as  Refs.  2  and  3  have  shown,  the  role  of  the  arc  is  simply  to  "plug  up"  a  fraction  of 
the  constrictor  exit  area,  with  most  of  the  flow  passing  around  the  arc.  The  temperature  of 
this  outer  gas  is  decoupled  from  the  arc,  and  is  governed  by  wall  heat  transfer,  particularly 
regenerative  heating  upstream  of  the  cathode,  and  this  temperature,  in  turn,  controls  the 
outer  flow  sonic  speed,  and  hence  the  specific  impulse.  The  detailed  iterations  arc  only 
now  being  conducted,  and  the  results  will  also  be  reported  in  the  upcoming  lEPC  paper. 
The  results  discussed  here  so  far  have  all  assumed  a  wall  temperature  of  lOOOK  upstream, 
increasing  to  1200K  at  the  constrictor  exit. 

j.  -NumierLcal  MQdsLfor  a  Hall  ThrusUr 

In  the  early  1960’s,  under  A. I.  Morozov,  the  Soviets  developed  the  first  concept  of 
a  thruster  with  closed  electron  drift  and  an  extended  acceleration  zone.  Between  1964- 
1969,  these  thrusters  were  improved  and  evolved  into  the  Hall  thruster.  Much 
experimental  work  has  been  performed  in  the  Soviet  Union  by  Morozov  (18,19,201^ 

Bugrova  I2i.22]_  Smirnov  1231^  Bishaev  124)^  Esipchuck  125]  and  Zubkov  126]  jq  list  a  few. 

A  great  deal  has  been  learned  from  the  data  about  thruster  design,  scaling,  plasma 
properties  in  the  acceleration  channel,  and  in  improving  performance. 
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Very  little  work  had  been  performed  on  numerically  modeling  these  thrusters. 
Kaufman  t^^l  modeled  performance  for  a  variety  of  assumed  potential  profiles.  Partial 
analysis  of  many  of  the  imponant  processes  were  reported  in  the  Soviet  literature  118-26) 
and  also,  more  recently  in  the  work  of  Yamagiwa  and  Kuriki  1-81.  Komurasaki  et.  al. 
at  the  University  of  Tokyo,  created  simplified  one  and  two  dimensional  models  for  the 
plasma.  These  models  assume  a  constant  ionization  to  loss  ratio,  use  a  simplified 
ionization  model,  and  neglect  thermal  conduction  in  the  electron  energy  equation. 

The  thrust  of  our  research  was  to  develop  a  more  detailed  and  accurate  model  of  the 
physical  processes  in  the  acceleration  region.  This  model  was  then  used  to  examine 
thruster  performance  over  various  operational  power  levels  using  Argon  and  Xenon  as 
propellants. 

This  research  uses  a  quasi-one  dimensional  transient  model  for  the  plasma 
properties  down  the  length  of  the  channel.  The  geometry  for  the  model  is  taken  from  an 
actual  thruster  for  comparison  of  the  numerical  results  with  collected  data.  The  magnetic 
field  is  assumed  to  be  imposed  and  taken  as  constant  in  time  and  the  plasma  is  assumed  to 
be  quasi-neutral.  Careful  treatment  is  given  to  the  electron  energy  equation,  including  heat 
conduction  in  its  derivation.  An  ion  distribution  function  as  a  function  of  the  electric  field, 
ionization  rates  and  wall  losses  is  defined  and  calculated.  Additionally,  momentum  transfer 
between  ions  and  neutrals  is  considered. 

The  leakage  of  electrons  from  the  downstream  cathode  to  the  upstream  anode  is  an 
imponant  process  in  the  physics  of  Hall  thrusters,  and  the  radial  magnetic  field  is 
introduced  to  control  this  leakage,  as  well  as  to  slow  down  the  flow  of  secondary  electrons 
to  the  anode.  This  increases  the  probability  of  a  given  electron  producing  ionizing  collision 
with  the  neutral  gas,  and  also  creates  a  high  density  electron  population  to  neutralize  the 
ions  during  their  acceleration.  We  have  modeled  the  plasma  impedance  to  the  flow  of 
electrons  by  an  extension  of  Bohm's  anomalous  diffusion  concepts,  following  Ref.  29. 
This  results  from  electron  scattering  by  turbulent  field  fluctuations.  The  Russian  literature 
121]  has  suggested  that  electron-wall  collisions  also  have  the  effect  of  interrupting  electron 
Larmor  gyrations,  thereby  inducing  additional  field-aligned  drifts  and  creating  a  near-wall 
shorting  layer.  This  effect  was  not  accounted  for  in  our  study.  The  results,  as  will  be 
shown,  appear  to  justify  our  approximations,  but  further  examination  of  the  possible 
effects  of  wall  shorting  is  needed. 
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Since  the  details  of  the  ni'^  ^iel  are  explained  in  the  attached  paper  AIAA  93-2491, 
we  will  here  only  comment  on  the  salient  aspects  of  the  results  and  on  perceived  remaining 
uncertainties  and  research  needs. 

Thw  internal  consistency  of  the  model  was  verified  by  careful  comparison  the  to  the 
data  of  Ref.  29  on  the  axial  distribution  of  (transversally  averaged)  ion  density,  electron 
temperature,  ionization  rate  and  potential  (see  Figs.  7-1 1  of  AIAA  93-2491),  which  are 
predicted  consistently  within  ±10%  of  the  data.  This  was  important  in  order  to  avoid 
compensating  errors  which  might  still  fortuitously  yield  good  performance  results. 
Especially  noteworthy  is  the  good  prediction  of  Te(x),  because  there  are  several  competing 
channels  for  the  energy  gained  by  electrons  from  the  field,  thermal  energy  being  only  one 
of  them.  This  Te,  in  turn,  determines  the  local  ionization  rate,  whose  accurate  prediction  is 
important  to  ensure  that  the  proper  efficiency  is  calculated. 

The  final  check  on  the  model  is  the  calculation  of  thruster  efficiency  when  only  the 
anode  current  and  the  mass  flow  rate  are  specified.  This  involves  calculation  of  the  thruster 
voltage  (not  specified)  and  the  beam  current  (smaller  than  the  anode  current  by  the  amount 
of  electron  upstream  leakage).  We  have  carried  out  computations  for  the  conditions  of  the 
experiments  of  Ref.  29,  using  Argon  propellant.  To  date,  we  have  only  completed  the  case 
of  m  =  lAeq,  over  the  tested  lanode  range.  The  results,  reproduced  in  Fig.  4  (from  our 
paper)  are  excellent.  Similar  calculations  covering 
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the  full  range  of  mass  flows  of  the  tests,  as  well  as  the  cases  w  here  Xe  was  used  as 
propellant,  are  now  underway,  and  will  be  reported  at  the  lEPC  Seattle  meeting  (Sept. 
1993).  Preliminary  results  are  very  encouraging. 

Despite  this  success,  several  areas  require  further  attention: 

(a)  As  noted,  wall  conductivity  was  not  included,  which  may  be  important  for  other 
conditions. 

(b)  The  loss  of  electron-ion  pairs  to  recombination  on  walls  was  modeled  originally 
by  assuming  perfect  plasma-wall  contact  and  using  the  flux  V  =  0.61n,i»,,  where 
Vg  is  the  Bohm  velocity.  This  proved  excessive,  and,  since  the  2-0  data  indicate 
that  the  plasma  near  the  anode  is  not  in  direct  contact  with  the  wall,  a  shape  factor 
was  included,  going  from  0  at  the  anode  to  1  at  mid-channel.  Clearly,  this  kind  of 
approximation  can  only  be  made  with  confidence  in  view  of  either  some  2-D  data 
or  some  2-D  computational  results. 

(c)  The  plasma  fluctuations  are  implicitly  accounted  for  in  the  anomalous 
resistance,  but  are  not  explicitly  calculated  otherwise,  nor  arc  any  other  potential 
macroscopic  instabilities  assessed.  Local  stability  analyses  were  reported  early  on 

but  these  can  only  serve  as  rough  guidelines,  and  much  more  work  is  needed 
here  (probably  in  terms  of  direct  simulations). 

(d)  Experimenters  report  wide  beam  divergence,  and  associated  lip  erosion  by 
impinging  ions.  This  is  another  2-D  effect  requiring  2-D  calculations. 

We  plan  to  work  to  these  problems  in  the  coming  two  years,  by  developing  a  2- 
dimensional  code  for  Hall  thrusters.  Initial  efforts  in  this  direction  are  aimed  at  assessing 
the  feasibility  of  directly  simulating  electron  transpon  by  means  of  particte-in-cell  codes, 
or,  alternatively,  extending  the  phenomenological  approach  used  so  far  to  the  2-D  case. 
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APPENDIX 

RELEVANT  PUBLICATIONS 


Following  are  the  papers  by  (a)  S.  Miller  and  M.  Martinez-Sanchez  on  arcjet 
fno« .  ling,  (b)  Ch.  Lentz  and  M.  Martinez-Sanchez  on  Hall  thruster  modeling,  (c)  E. 
Niewood  and  M.  Martinez-Sanchez  on  anode  drops  in  MPD  thruster,  and  (d)  M.  Martinez- 
Sanchez  and  A.  Sakamoto  on  a  simplified  arcjet  model.  These  contain  detailed  technical 
accounts  of  the  topics  discussed  above. 
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Abstract 

A  detailed  model  has  been  developed  to  study 
the  gasdynamic  flow  in  an  electrothermal  arcjet 
thruster  through  numerical  simulation.  This  two- 
temperature  model  consistently  incorporates  viscos¬ 
ity  and  heat  conduction  effects  by  expressing  the 
governing  equations  in  a  Navier-Stokes  form.  Other 
simulated  effects  include  ohmic  dissipation,  coUi- 
sional  energy  transfer  between  electrons  and  heavy 
species,  ambipolar  diffusion,  nonequilibrium  dissoci¬ 
ation  and  ionization,  and  continuum  radiation.  The 
fluid  equations  are  solved  by  MacCotmack’s  method, 
while  an  iterative  procedure  is  used  to  relax  an  elec¬ 
tric  potential  equation,  from  which  the  current  dis¬ 
tribution  in  the  thruster  is  obtained.  Converged  so¬ 
lutions  are  compared  with  experimental  results  from 
the  German  TTl  radiatively-cooled  arcjet  thruster 
with  hydrogen  propellant.  Results  are  presented  for 
a  baseline  case  which  reveal  the  two-dimensional, 
two-fluid  nature  of  the  interior  flow,  especially  in 
terms  of  the  distribution  and  anode  attac*  m.nt  of 
the  electric  current  and  the  growth  and  development 
of  the  arc  region.  Calculated  discharge  voltage  is 
within  a  few  percent  of  experimental  measurements, 
and  predicted  specific  impube  b  within  5-10%  agree¬ 
ment  over  a  range  of  operating  parameters. 

Nomenclature 

c.  Specific  heat  at  constant  volume  [J [moUI^K] 
Da  Ambipolar  diffusion  coefficient  [m^/s] 

d, r  Ambipolar  flux  of  ions  and  electrons  [l/m^/«] 

e  Electric  charge  [C],  internal  energy  [J/kg\ 

E  Electric  field  [V/m] 

Ei  Ionization  energy  [J] 

Ei  Dbsociation  energy  [J] 

El  Elastic  collbional  energy  transfer  [W/m^] 

Eaii  Vibrational  excitation  energy  [J] 
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h  Planck’s  constant  [Js],  enthalpy  [J /kg] 

I  Total  current  [A] 

j  Electric  current  density  [A/m*] 

k  Boltsman’s  constant  [J/^K] 

K  Equilibrium  constant 

m  Particle  mass  [kp] 

n  Number  density  [l/m®] 

li.  Net  electron  production  rate  {l/m*/s] 

nj  Net  atomic  hydrogen  production  rate  [l/m*/s] 

N  Avogadro’s  number  [1/moie] 

p  Scalar  pressure  [Pa] 

q  Heat  flux  vector  [TV/m*] 

B.  Real  gas  constant  [Jlkgl^K] 

R  Universal  gas  constant  [J/mole/*K] 

R  Energy  loss  due  to  coatinaum  radiation  [W/rr?] 

T  Temperature  [*JC] 

u  Mean  flow  velocity  [m/s] 

vs  Bohm  velocity  [m/s] 

V  Slip  velocity  [m/s] 

X  Mole  fraction 

a  Ionization  fraction 

K  Coefficient  of  thermal  conductivity  [W/m/*K] 

ft,  Coefficient  of  viscosity  [kp/ms] 

u  CoUbion  frequency  [1/s] 

4  Electric  potential  [V] 

^  Viscous  dissipation  function  [W/m*] 

V*  Electron  mobility  [m*/n/C] 

p  Mass  density  [kp/m®] 

<T  Electrical  conductivity [m/io/m],  molar  concentration 
Q  Average  effective  coUbion  integral  [m*] 


1  Introduction 

Low  power  arcjet  thrusters  have  recently  been  flight 
qualified  through  ground  testing  and  will  soon  be 
tested  in  space  for  stationkeeping  of  geosynchronous 
communication  satellites.  Most  of  the  impetus  for 
design  strategies,  however,  has  come  ftom  empir¬ 
ical  studies  and  experimentation,  and  a  need  re¬ 
mains  to  better  understand  the  underlying  physics, 
detailed  energy  balances,  and  transport  mechanbms 
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of  these  devices.  In  ordet  to  fully  understand  the 
physics  of  aicjet  thrusters,  it  is  necessary  to  accu¬ 
rately  model  many  complex  interacting  phenomena. 
Unfortunately,  the  complexity  of  the  modeb  and 
equations  needed  to  accurately  represent  these  de¬ 
vices  effectively  limits  the  use  of  analytic  techniques 
to  simplified  cases,  through  which  one  may  obtain 
useful  physical  insights  but  inadequate  predictions 
of  thruster  performance.  Experimental  techniques 
provide  much  useful  empirical  data,  but  many  quan¬ 
tities  of  interest  are  not  accessible  in  the  impor¬ 
tant  regions  of  the  thrusters.  For  these  reasons  nu¬ 
merical  methods  of  solving  the  governing  equations 
have  become  an  important  tool  for  conducting  arejet 
thruster  research. 

Previous  numerical  modeling  of  arejet  thrusters 
has  focused  on  the  development  of  1-D,  2-D,  and 
axbymmetric  modeb  with  relatively  simple  physics 
and  geometries.  The  level  of  detail  has  ranged 
from  1-D  modeb[l]  to  coupled  quasi-analytic  mod¬ 
eb  of  the  inner  (arc)  and  outer  (cold  gas)  flows[2, 
3,  4],  to  simplified  axbymmetric  space-marching 
techniques[5],  and  finally  to  2-D  and  axbymmet¬ 
ric  viscous  codes  which  begin  to  incorporate  most 
of  the  detailed  physical  processes[6,  7,  8].  The  lat¬ 
est  research  has  obtained  results  which  variously  in¬ 
clude  ohmic  heating,  electron  heat  conduction,  ther¬ 
mal  and  ionbational  nonequilibrinm,  and  empbical 
modeb  of  radiation  losses.  There  are  still  a  num¬ 
ber  of  bsues,  such  as  vbcous  and  diffusive  effects, 
arc  formation  and  attachment,  the  heat  balance  in 
the  anode,  and  ultimately  the  accurate  prediction  of 
voltage  and  efficiency,  which  need  to  be  addressed. 
Thb  paper  describes  a  generalbed,  more  physically 
accurate  model  of  the  gasdynamic  flow  through  an 
arejet  thruster,  which  includes  the  aforementioned 
effects  and  compares  favorably  to  experimental  re¬ 
sults  obtained  with  medium  power  hydrogen  arejets. 

2  Model 

2.1  Basic  Assumptions 

A  basic  diagram  of  an  electrothermal  arejet  thruster 
b  shown  in  Figure  1.  Thb  model  b  based  on  an 
aixbymmetric  formulation.  Consequently,  variations 
in  flow  quantities  in  the  azimuthal  direction  are  ne¬ 
glected.  A  component  of  the  flow  velocity  in  the 
^-direction,  however,  b  incorporated  to  account  for 
the  “swirl”  injection  of  most  experimental  arejets. 
Thb  injected  vortex  has  been  shown  to  stabilize  arc 
attachment  and  to  help  maintain  a  concentrated  arc 
along  the  centerline  of  the  thruster.  Although  u«  b 
assumed  constant  in  the  fl-direction,  it  is  allowed  to 
develop  in  the  axial  and  radial  directions  through 


the  conservation  of  global  angular  momentum. 

The  model  has  been  developed  in  a  general  enough 
sense  so  that  any  monatomic  or  diatomic  propellant 
may  be  simulated.  For  the  purpose  of  thb  research 
hydrogen  was  selected  as  the  propellant  of  choice, 
due  to  its  low  molecular  weight  (and  therefore  high 
performance)  and  its  simple  molecular  structure, 
which  allows  for  analytic  evaluation  of  the  necessary 
transport  coefficients.  Nonequilibrium  dbsociation 
and  ionization  are  modeled,  and  four  species  of  par¬ 
ticles  are  tracked:  diatomic  molecules,  monatomic 
neutrab  and  ions,  and  electrons.  Dissociation  b 
modeled  by  heavy  species  coUbions  and  by  electron 
impact,  and  the  ionbation  process  b  based  on  elec¬ 
tron  impact  ionbation  and  three-body  recombina¬ 
tion,  with  only  H*  ions  considered. 

The  foUowing  assumptions  are  made  regarding  the 
state  of  the  flow  in  the  thruster  and  the  physical  pro¬ 
cesses  involved.  The  plasma  produced  by  ionbing 
electron  coUbions  b  assumed  to  be  macroscopically 
neutral,  so  that  n«  =  n^.  Strong  coupling  b  assumed 
between  the  ions  and  neutrab,  designated  together 
as  the  heavy  species.  Thb  implies  that  u,-  u),  2:  u 

(except  for  ambipolai  diffusion),  and  Tj  2  Tn  S  Tg. 
Effects  which  axe  consbtently  incorporated  include 
ambipolat  diffusion,  heat  conduction,  viscous 
and  dissipation,  ohmic  heating,  coUbionai  energy 
transfer  between  electrom  and  heavy  species,  and 
energy  lost  through  continuum  radbtion.  The  self- 
induced  magnetic  field  of  the  ionbed  gas  b  neglected 
due  to  the  low  current  density  and  thus  negli^ble 
magnetic  pressure  in  arejet  thrusters,  and  the  indi¬ 
vidual  species  are  assumed  to  obey  the  ideal  gas  law. 
Given  the  aforementioned  assumptions,  the  model 
can  be  summarbed  by  a  set  of  nine  partial  differen¬ 
tial  equations  which  must  be  solved  locally  in  order 
to  generate  a  viable  simulation  of  the  flow  in  an  ar¬ 
ejet  thruster. 


The  set  of  equations  which  govern  the  flow  in  the 
model  arejet  thruster  of  thb  research  b  essentially 
a  group  of  modified  Navier-Stokes  equations.  These 
include  an  equation  of  state  and  equations  for  the 
ion,  neutral  atom,  and  global  density;  the  axial,  ra¬ 
dial,  and  azimuthal  global  momentum;  the  electron 
and  heavy  species  energy;  and  the  electric  potential. 
In  order  to  maintain  numerical  robustness,  these 
equations  are  written  in  as  conservative  a  manner 
as  possible,  particularly  with  respect  to  the  i  ternu 
which  appear  due  to  the  axbymmetry  of  the  prob¬ 
lem. 

The  electric  potential  equation  b  derived  by  com¬ 
bining  Ohm’s  law,  j  =  with  the  equation 


V  ■  j  =  0,  where  ^  = 


is  the  electron  mo¬ 


bility.  Assuming  a  potential  of  the  form  E  =  -V0, 
the  resulting  equation  is  given  by 

d  ,  ,  /d*p.  i5p,  a*p.\ 


where 


dp,  av»  ap,  aV’ 

^  dr  dr  dz  dz  ' 


,  ,  H 

and  4),  =  <r  —  . 


The  current  density  may  then  be  extracted  by  solv¬ 
ing 

dp,  d4>  ,  .  dp,  d<t> 

The  global  density  equation  is  obtained  by  sum¬ 
ming  the  individual  species  continuity  equations: 

dpr  dpu,r  d(ni,r 

IT  +  “sT  +  “sT  =  “■  <’> 

For  the  ion  (or  electron)  density,  the  governing  equa¬ 
tion  is  derived  from  the  statement  of  ion  mass  con¬ 
servation,  modified  to  account  for  the  ambipolar  flux 
of  charged  particles  in  terms  of  ion  density  gradients: 

dp,r  d{p,Ur+d,r)r  d{p,u,  +  d,,)r 

-W  * - ^ - Vr - = 

(5) 

where  d,  is  the  ambipolar  flux,  given  by 

V  ^TTlyg  Qiny^t  ^n) 

(6) 

The  source  term  n,  represents  the  net  rate  of  pro¬ 
duction  of  ions  per  unit  volume  through  inelastic 
collisional  processes.  The  statement  of  neutral  atom 
mass  conservation  is  obtained  in  a  similar  man¬ 
ner,  and  the  effect  of  ambipolar  diffusion  is  con¬ 
sistently  included  by  assuming  that  both  neutral 
species  travel  at  the  same  velocity  and  by  employ¬ 
ing  the  definition  of  the  mean  flow  velocity  (u  = 

dpH^  dpH^rr  g  /  PH  \ 

dt  ^  dr  dr  \(pn,  +  pa)  m,  ) 

^  dpHU,r  d  /  Pa 

dz  dz  V(PJ?,  +Ph)  m,  “  ) 
^ma(na-*r  <<rv>n,na.,-n,)r.  (7) 

Here  the  source  term  represents  the  net  rate  of  pro¬ 
duction  of  neutral  atoms,  given  by  the  difference  be¬ 
tween  the  net  rate  of  production  due  to  dissociation 


and  the  net  rate  of  ionization  of  those  neutral  atoms 
produced. 

The  three  global  momentum  equations  are  given 
by 

dpu,r  d(pnl  +p-r,,  )r  ^(pu^u,  -  r„)r 
di  dr  dz 

-I-  r„  -  pui  -  p  =  0,  (8) 

dputr  d(pu,u#  -  r,0)r  5(pusu,  -  r,,)r 

dt  dr  dz 

+  pUrUt  =  0,  and  (9) 

dpU,r  d{pUrU,  ,  ^(P“? +P- tr,)r 

+  _  +  - -  =  0; 

(10) 

and  the  corresponding  viscous  stress  tensor  compo¬ 
nents  in  axisymmetric  coordinates  are 


dUr 

dUg 

(11) 

TVf 

=  (2 

dr 

dz 

■t) 

2 

dUr 

dug\ 

(12) 

T$$ 

= 

'  r 

dr 

2  Yo 

du. 

dUr 

(13) 

T„ 

=  3^^' 

~d^  ■ 

■  17  ■ 

■t) 

(du. 

V 

(14) 

TV#  =  P, 

[dr 

T*  > 

) 

(dUr 

dUg 

(15) 

TV*  =-  Mf 

■''17 

) 

dtta 

(16) 

T»*  : 

II 

?|, 

The  heavy  species  energy  equation  is 

^Pgs^  g(pg>ur  +  g>r)v  9(pefU,+q,,)r 
dt  dr  dz 

^[PH,ea,Va,r  +  PgegVg,  +  Pa+ea+Va^,)r 
dr 

=  (* + 

(17) 

where  the  internal  energy  and  enthalpy  are  defined 
by 
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pe,  =  -  (pH  +  Ph+  )  RhT,  +  -^PaiRBiT, 

+  -  Lr«t,  CO 

e^  —  1  ^ 


phf  =  peg  +  pa,RH,Tj  +  (pa  +  Pa*)  RaT,.  (19) 
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Additional  flux  terms  have  been  included  in  the  ra¬ 
dial  direction  to  account  for  the  radial  transport 
of  energy  by  ambipolar  diffusion;  these  terms  are 
dependent  on  the  internal  energy  and  slip  veloc¬ 
ity  (V^,)  of  each  species.  The  heat  flux  vectors  in 
Eqn.  17  can  be  expressed  as 

9*.  = 

while  the  viscous  dissipation  function  is  given  by 


The  collisional  transfer  of  energy  from  electrons  to 
the  heavy  species  is  represented  by  the  following  ex¬ 
pression: 

El  ~  3-^(i/eir+  +  u,B  +  “  'Ef)'  (22) 

ntH 

The  coefficient  6,  in  Eqn.  22  is  necessary  to  correct 
for  the  fact  that  electron- collisions  are  inelastic  in 
nature.  An  internal  energy  form  of  the  heavy  species 
energy  equation  is  utilised  rather  than  a  more  con¬ 
servative  total  energy  form  because  of  numerical  con¬ 
cerns.  In  some  regions  of  an  arcjet  thruster  the  ki¬ 
netic  energy  and  the  net  loss  of  energy  due  to  dissoci¬ 
ation  can  become  so  large  that  numerical  difficulties 
arise  if  these  terms  are  kept  in  a  total  energy  state 
vector.  The  same  is  true  of  the  net  loss  of  energy 
due  to  ionisation  in  the  electron  energy  equation. 

The  conservation  of  energy  for  electrons  is  given 
by  the  following  expression: 

dp,  E,t  d{p,  u,r  H,  +q,r)r  d  (p,  u,,  ff,  +q„)r 

dt  dr  dz 


The  final  equation  requited  for  closure  of  the  set 
is  the  equation  of  state: 

P=  Sp>  =  +  (26) 
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2.3  Dissociation  and  Ionization 

The  nonequilibrium  dissociation  rate  hjf  is  derived 
by  following  the  procedure  and  nomenclature  of 
Biasca[9]  for  collisions  of  heavy  species.  Accordingly, 
this  rate  is  given  by 

hg  =  ANT^  exp  {.^a<rs  +  mg^ffg,) 

A 

"  f 

where  N  is  Avogadro’s  number,  is  the  equilib¬ 
rium  constant  in  terms  of  partial  pressures,  the  <r^s 
are  the  species  molar  concentrations,  and  the  ap¬ 
propriate  constants  for  hydrogen  are  listed  in  Ta¬ 
ble  1(10].  Dissociation  by  e  -  collisions  is  repre¬ 
sented  by  the  second  term  on  the  right-hand  tide  of 
Eqn.  7,  where  the  reaction  rate  coefficient  <  av  > 
as  a  function  of  T*  is  taken  from  Janev  et  al.[ll]. 

Table  1:  Constants  for  the  Hydrogen  Dissociation 
Rate 


For  ionization,  the  finite  production  rate  is 
given  by  the  generalised  model  of  ionisation 
and  three-body  recombination(12]  as  modified  by 
Sheppard[13]: 


=  - El  -  Ei  <  <Tv  >  n,ng,  -  Eih,  -  RJ  r. 

(23) 

Here  ^  represents  heating  due  to  ohmic  dissipation, 
R  is  the  radiative  loss,  and  the  total  energy  and 
enthalpy  are  given  by 

E,  =  ^R.T.  +  (24) 

and 

H.  =  E,  +  R,T,.  (25) 

The  required  electron  velocities  are  extracted  from 
the  local  current  density. 


n,  =  Rn,  {Sng  -  nj) 


(») 


«  =  6.985  X  10““e*p 


2.4  Transport  Properties 

Because  of  the  multi-component  nature  of  the  gas 
in  an  electrothermal  arcjet,  even  with  hydrogen  ax 
a  propellant,  the  equations  for  the  transport  prop¬ 
erties  become  quite  complex.  Since  data  are  more 
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readily  available  for  hydrogen  in  the  form  of  col¬ 
lision  integrals,  a  formulation  of  the  transport  co¬ 
efficients  based  on  these  integrals  rather  than  on 
collision  cross-sections  is  implemented.  The  kinetic 
theory  of  multicomponent  gases  has  been  rigorously 
applied  to  the  calculation  of  transport  coefficients 
by  Curtiss  and  Hirschfelder[14],  and  the  formulation 
here  follows  their  work. 

For  the  electrical  conductivity,  a  first  order  ap¬ 
proximation  by  Grier[l5]  to  the  rigorous  kinetic  the¬ 
ory  has  been  chosen  for  use  in  this  research: 


°  16  V  irrujAiT, 


^(1,1)  ’ 

where  the  Xj  are  mole  fractions  and  the  are  av¬ 
erage  effective  collision  integrals.  The  variable  At  is 
a  more  complicated  function  of  mole  fractions  and 
additional  collision  integrals. 

The  gas  coefficient  of  viscosity  is  calculated  based 
on  a  mean  free  path  mixture  formula,  which  is  a 
function  of  the  collision  integrals  and  the  species 
number  densities  and  pure  viscosities[12]; 
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+'‘HV3^ 


na  +  ng 
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i  +  2«h+ 

'‘h-H 


+  2ng 

where  the  pure  viscosities  are  given  to  lowest  order 
by 

=  2.6693  X  (34) 

Electron  and  heavy  species  thermal  conductivity  co¬ 
efficients  are  calculated  based  on  similar  mean  free 
path  arguments  and  mixture  rules,  and  the  pure 
thermal  conductivities  are  given  by  Euken’s  relation 

15  /4  c,  3^ 

where  c,  equals  |  R  for  monatomic  gases  and  |  R  for 
diatomic  gases. 

Collision  integrals  required  in  the  calculation  of 
transport  coefficients  are  interpolated  as  a  function 


of  temperature  from  data  by  Belov[l6j,  Grier[15], 
and  Vanderslice  et  al.[17].  The  accuracy  of  the  above 
approximate  formulas  for  the  transport  coefficients 
was  verified  by  comparison  to  previous,  more  de¬ 
tailed  work  for  hydrogen  in  thermal  and  ionizational 
equilibrium. 


3  Numerical  Method 

3.1  Integration  Scheme 

The  numerical  method  used  in  this  research  has  been 
previously  implemented  in  numerical  simulations  of 
MPD  thrusters  with  coaxial  geometries[18].  The 
governing  fluid  equations  are  numerically  integrated 
using  MacCormack’s  explicit  node-based  method. 
This  predictor-corrector  method  is  second  order  ac¬ 
curate  in  both  time  and  space,  and  it  has  proved  an 
excellent  means  of  solving  the  compressible  Navier- 
Stokes  equations  for  a  variety  of  conditions.  The 
explicit  nature  of  the  algorithm  allows  for  simpler 
coding  than  an  implicit  scheme  because  of  the  com¬ 
plex  nature  of  the  source  terms  and  boundary  con¬ 
ditions,  and  it  also  provides  a  means  of  obtaining 
time-accurate  solutions. 

An  empirical  stability  formula  developed  by  Mac- 
Cormack  and  Baldwin[19]  is  utilised  to  calculate  the 
integration  time  step,  which  is  then  mnltiplied  by 
a  fractional  coefficient  to  account  for  the  effect  of 
source  terms  and  nonlinearities.  The  time  accuracy 
of  the  solution  is  maintained  so  long  as  the  same  time 
step  is  used  at  all  computational  grid  locations.  Al¬ 
though  MacCormack’s  method  contains  some  inher¬ 
ent  dissipation,  additional  numerical  smoothing  is 
required  in  order  to  damp  unwanted  numericsd  oscil¬ 
lations  in  the  fluid  equations.  Consequently,  simple 
2nd  and/or  4th  order  smoothing  terms  are  applied 
to  each  equation  as  necessary,  based  on  the  method 
of  Kutler,  Sakeli,  and  Aidlo(20]. 

The  electric  potential  equation,  being  predomi¬ 
nantly  elliptic  in  nature,  is  solved  by  iteration  nring 
a  successive  overrelaxation  (SOR)  technique.  Since 
the  physical  grid  is  nonuniform  in  order  to  closely 
represent  the  geometry  of  the  actual  arejet  being 
modeled,  the  governing  equations  are  transformed 
into  natural  coordinates  and  then  solved  on  a  uni¬ 
form  Cartesian  computational  mesh.  The  physical 
grid  used  for  the  simulations  presented  in  this  paper 
is  shown  in  Figure  2. 

3.2  Boundary  Conditions 

The  conditions  at  the  inlet  of  the  computational 
domain  are  postulated  to  be  essentially  those  of  a 
flow  which  has  just  been  injected  into  the  thruster 
plenum  by  a  large  number  of  evenly  spaced  jets.  The 


flow  is  therefore  assumed  to  be  subsonic  and  paral¬ 
lel  to  the  thruster  walls.  A  Gtaction  of  the  total  inlet 
velocity,  typically  30-50%,  is  specified  as  being  in 
the  asimuthal  direction  in  order  to  simulate  an  in¬ 
jected  swirl.  The  mass  flow  rate  and  total  enthalpy 
are  specified  based  on  the  particular  run  parame¬ 
ters,  and  the  ionization  fraction  is  set  to  a  small 
value,  typically  1  x  10"®.  The  density  is  obtained 
from  a  downwind  finite  difference  approximation  of 
the  overall  continuity  equation,  and  the  axial  veloc¬ 
ity  i.s  then  found  from  the  specified  mass  flow  rate. 
The  inlet  electron  temperature  is  set  equal  to  that 
of  the  next  inside  point,  and  the  heavy  species  tem¬ 
perature  is  then  found  from  the  total  mthalpy.  No 
current  b  allowed  to  leave  the  domain  at  the  inlet 
(J.  =  0). 

The  boundary  conditions  at  the  outlet  of  the 
thruster  depend  on  whether  the  exit  flow  is  subsonic 
or  supersonic.  In  both  cases  the  electron  tempera¬ 
ture  is  set  equal  to  that  of  the  next  inside  point,  and 
no  current  may  pass  beyond  the  exit  plane.  If  the 
flow  is  supersonic  at  a  point  on  the  exit  plane,  then 
the  remaining  quantities  are  extrapolated  from  their 
values  at  the  preceding  two  grid  points  of  the  mesh. 
If  the  flow  is  subsonic,  then  the  exit  pressure  is  set 
equal  to  a  small  value  representing  neat- vacuum  con¬ 
ditions  and  the  density  and  axial  velocity  ate  given 
by  Riemann  invariamts.  The  remaining  quantities 
at  a  subsonic  outlet  are  then  calculated  as  in  the 
supersonic  case. 

For  those  boundary  points  lying  beyond  the  tip  of 
the  cathode  on  the  line  of  symmetry  (r  =  0),  the 
radial  and  azimuthal  flow  velocities  ate  set  equal  to 
zero  and  a  zero  radial  gradient  is  imposed  on  the 
remaining  quantities. 

Flow  boundary  conditions  at  the  thruster  wails 
include  viscous  no-slip  conditions  for  the  axial  and 
radial  fluid  velocities.  The  electron  temperature  is 
set  equal  to  that  of  the  next  inside  grid  point,  and 
the  heavy  species  temperature  is  held  constant  at 
1000°K  upstream  of  the  constrictor,  increasing  lin¬ 
early  to  nOO^K  at  the  constrictor  exit.  This  profile 
was  chosen  based  on  experimental  and  numerical  cal¬ 
culations  of  the  anode  wall  temperature  distribution 
for  a  reasonable  operating  range  of  the  German  TTl 
radiatively-cooled  arcjet  thruster[21].  On  the  cath¬ 
ode  the  wall  temperature  is  allowed  to  increase  to 
a  maximum  of  2000*K  at  the  tip.  For  the  bound¬ 
ary  condition  on  electron  density  at  each  electrode, 
a  balance  is  postulated  between  the  flux  of  ions  ar¬ 
riving  at  the  sheath  edge  by  ambipolar  diffusion  and 
the  flux  of  ions  arriving  at  the  wall  by  virtue  of  their 
thermal  energy  at  the  Bohm  velocity  (t)B)[22]: 

=  0.37n,Wfl,  (36) 

dy 


where 


V  m. 


(37) 


This  boundary  condition  neglects  the  voltage  drops 
present  in  the  non-neutral  plasma  sheath.  Also 
omitted  is  the  dependence  of  the  thermionic  emis¬ 
sion  of  electrons  at  the  cathode  on  the  cathode  wall 
temperature.  Application  of  the  perpendicular  over¬ 
all  momentum  equation  at  the  walls  provides  the 
approximate  condition 


dp 

^l«.H=0,  (38) 

where  inertial  and  viscous  terms  have  been  ne¬ 
glected. 

The  boundary  condition  on  the  wail  electric  po¬ 
tential  is  that  there  is  no  current  perpendicular  to  an 
insulating  section,  and  the  potential  on  the  anode  is 
set  equal  to  a  fixed  but  arbitrary  voltage.  For  numer¬ 
ical  reasons,  anode  current  attachment  is  restricted 
to  that  portion  of  the  outer  wall  downstream  of  the 
constrictor  exit.  On  the  cathode  tip,  a  uniform  ax¬ 
ial  current  density  is  prescribed  which  sums  to  the 
specified  total  current;  the  potential  at  the  cathode 
is  then  chosen  so  as  to  maintain  this  current  level. 
A  cathode  voltage  drop  equal  to  the  ionisation  po¬ 
tential  plus  one  half  of  the  dissociation  potential  of 
the  gas  is  added  to  the  calculated  voltage  in  order 
to  account  for  the  model’s  neglect  of  the  cathode 
sheath  region.  In  addition,  an  anode  voltage  drop 
is  subtracted  from  the  calculated  voltage  in  order  to 
account  for  the  anode  sheath  region.  A  negative  po¬ 
tential  gradient  is  required  in  this  sheath  in  order  to 
turn  back  excess  electrons  since  the  extracted  cur¬ 
rent  is  much  less  than  the  random  thermal  flux  of 
electrons  to  the  anode  wall  ^  ^cn,c,).  This 

voltage  drop  is  given  by 


3.3  Procedure 

Once  the  initial  conditions  of  the  flow  have  been 
specified,  the  solution  is  calculated  numerically  as 
follows.  First,  the  timesteps  are  calculated  for  the 
overall  and  heavy  species  equations  at  ail  grid  points 
from  stability  criteria.  The  most  restrictive  time 
step  is  chosen  as  the  time  interval  for  the  current 
integration  step,  and  these  equations  are  integrated 


according  to  MacCormack’s  method.  Next,  the  elec¬ 
tron  density  and  energy  equations  are  integrated 
based  on  their  individual  stability  restrictions  while 
maintaining  consistency  with  the  previously  deter¬ 
mined  time  interval.  The  electron  equations  are  in¬ 
tegrated  separately  because  the  associated  stability 
criteria  can  be  up  to  30  times  more  restrictive  than 
for  the  overall  flow  equations.  Finally,  the  electric 
potential  equation  is  relaxed  a  specified  number  of 
iterations.  This  system  of  equations  is  then  repeat¬ 
edly  integrated  until  the  solution  is  considered  con¬ 
verged.  Computations  were  performed  on  a  vari¬ 
ety  of  machines,  including  DECStation  5000s,  IBM 
RS-eOOOs,  Silicon  Graphics  bis  bidigos,  and  Cray 
X-MPs. 


4  Results 

Results  have  been  achieved  for  comparison  to  the 
German  TTl  radiation-cooled  arcjet  thruster  at  a 
number  of  operating  points.  For  the  baseline  case 
of  /  =  1004  and  m  =  O.lj/s,  Table  2  compares 
bulk  results  of  the  simulation  to  experimental  mea¬ 
surements  and  Figures  3  through  10  show  line  and 
contour  plots  of  representative  quantities.  There 
is  excellent  agreement  in  the  voltage  between  pre¬ 
dicted  and  experimental  results,  and  good  agree¬ 
ment  (within  7%)  in  the  thrust  and  specific  impulse. 
The  accuracy  in  voltage  prediction  shows  that  the 
model  is  accurately  modeling  arc  growth,  electrical 
conductivity,  and  current  attachment  given  the  as¬ 
sumptions  we  have  made  regarding  electrode  sheath 
voltage  drops.  A  major  reason  for  the  discrepancy 
in  thrust  prediction  b  that  the  specified  inlet  gas 
temperature  is  probably  too  high.  Recent  modeling 
of  the  anode  heat  balance  in  this  thruster  has  shown 
that  the  inlet  temperature  for  the  baseline  case  b 
about  875‘’K  rather  than  1000*if  as  specified  in  the 
present  boundary  conditions.  Since  the  thrust  scales 
as  the  plenum  total  pressure,  which  scales  as  y/7o, 
a  12%  decrease  in  the  inlet  gas  temperature  could 
translate  into  a  6%  reduction  in  thrust.  Thb  would 
bring  the  simulation  results  into  very  close  agree¬ 
ment  with  the  experimental  data. 


Table  2:  Comparbon  of  Predicted  and  Experimental 

Results  for  Baseline  Case - - 

p-  I  Predicted  ]  Experiment  | 


Voltage  (V) 

114 

112  ! 

Power  (kW) 

11.4 

11.2 

Thrust  (N) 

1.01 

0.94 

Specific  Impube  (s) 

1030 

960 

Efficiency 

0.442 

0.395 

Figure  3  shows  an  axial  line  plot  of  the  electric 
potential  bom  the  cathode  tip  to  the  anode  attach¬ 
ment  zone.  The  neax-cathode  voltage  drop  AK,  b 
composed  of  a  15.8V  drop  assigned  to  the  cathode 
sheath  (Vi  -t-  jVj)  and  an  8V  drop  calculated  in  the 
first  few  gridpoints  downstream  of  the  tip.  The  near¬ 
anode  voltage  drop  AV.  %  15V  b  composed  of  a 
22V  drop  captured  by  the  simulation  and  a  -7V  drop 
associated  with  the  electron-repelling  sheath.  Thb 
total  anode  voltage  drop  b  associated  with  the  net 
deposition  of  energy  into  the  anode  block  by  heavy 
species  heat  conduction  and  by  the  impingement  of 
current-carrying  electrons.  Assuming  that  the  en¬ 
ergy  transferred  per  unit  area  b  of  the  form 

dT  r  5  1 

2^.  +  ' 

(40) 

using  the  results  of  the  baseline  flow  simulation 
yields  an  equivalent  voltage  of  14.5V  for  thb  de¬ 
posited  power,  which  agrees  well  with  the  AV,  seen 
in  the  potential  profile. 

Current  streamlines  are  plotted  in  Figure  4.  In 
thb  case,  the  bulk  of  the  current  attaches  within  the 
first  quarter  of  the  notsle,  with  a  peak  just  down¬ 
stream  of  the  constrictor  exit.  The  flow  becomes 
fuUy  ionued  along  the  centerline  immediately  down¬ 
stream  of  the  cathode  tip  and  remains  so  through 
the  first  part  of  the  nossle  expansion,  beyond  which 
there  b  some  recombination  (Figure  5).  The  bound¬ 
ary  of  the  partially  ionued  region  grows  to  approx¬ 
imately  50%  of  the  channel  by  the  constrictor  exit, 
and  thb  region  b  essentially  entrained  in  the  flow 
throughout  the  nossle.  The  primary  heating  mecha- 
nbm  b  ohmic  dbsipation,  which  peaks  locally  along 
the  constrictor  centerline  and  just  beyond  the  con¬ 
strictor  exit  near  the  anode.  Thb  b  evidenced  by 
the  local  maxima  in  electron  temperature  in  these 
regions  (Figure  6). 

Within  the  highly  ionued  region  of  the  arc  in  the 
constrictor,  collbional  energy  transfer  between  elec¬ 
trons  and  heavy  species  rabes  the  gas  temperature  to 
20, 000  -  30,  OOO'/ir,  or  neatly  the  same  temperature 
as  the  electrons.  Outside  of  the  arc  the  flow  remains 
cold,  at  a  temperature  approximately  equal  to  the 
anode  wall  temperature.  Thb  outer  flow  remains 
essentially  uncoupled  from  the  hot  core  flow,  as  pic¬ 
tured  in  Figure  7.  The  flow  velocity  b  also  uncou¬ 
pled,  although  vbcous  forces  eventually  wash  out  the 
separation  in  the  nozzle  expansion  (Figure  8).  Since 
the  pressure  b  nearly  uniform  in  the  radial  duection 
(Figure  9),  rapid  acceleration  of  the  core  flow  occurs 
throughout  the  low  pressure  region  of  the  arc  in  the 
constrictor.  Once  the  bulk  of  the  pressure  work  has 
been  utilized  in  the  expansion  process,  however,  vis¬ 
cous  forces  arbing  from  steep  velocity  gradients  in 


Table  3:  Average  Ratio  of  Terms  to  Oominant  Term 
in  Electron  Energy  Equation  _ 


TERM 

ARC  1  OUTER  1 

Radial  Convection 

0.338 

0.423 

Axial  Convection 

0.528 

0.444 

Radial  Conduction 

0.400 

0.216 

Axial  Conduction 

0.007 

0.056 

Ohmic  Dissipation 

0.508 

0.504 

Transfer  to  Heavy  Species 

0.023 

0.069 

Net  Dissociation 

0.006 

0.007 

Net  Ionization 

0.345 

0.207 

Radiation 

0.001 

0.001 

Table  4;  Average  Ratio  of  Terms  to  Dominant  Term 
in  Electron  Density  Equation  _ 


1  TERM 

ARC  I 

OUTER 

;  Radial  Divergence  I 

0.604 

0.636 

Axial  Divergence 

0.830  1 

0.788 

1  Radial  Ambipoiar  Diffusion 

0.512  1 

0.212 

.  -  -j 

!  Axial  Ambipoiar  Diffusion 

0.019  i 

0.009 

Net  Ionization 

1  0.240  i 

0.301 

the  central  core  decelerate  the  flow  signiflcantly  in 
the  nozzle  expansion.  Both  the  inner  and  outer  flows 
accelerate  smoothly  through  sonic  velocity  at  and 
just  beyond  the  constrictor  exit.  Figure  10  shows 
Mach  number  contours  for  the  baseline  case  arcjet 
simulation.  Interestingly,  due  to  the  large  radial  gra¬ 
dients  in  temperature,  composition,  and  velocity,  the 
peak  Mach  number  occurs  in  the  interior  region  of 
the  flow  rather  than  at  the  centerline. 

Simulating  the  current  attachment  at  the  anode 
realistically  and  self-consistently  has  been  a  major 
difficulty  in  previous  arcjet  simulations.  The  effec¬ 
tiveness  of  this  model  in  simulating  this  region  is 
due  to  the  incorporation  of  separate  energy  equa¬ 
tions  for  the  heavy  species  and  electrons  and  to  the 
use  of  nonequilibrium  dissociation  and  ionization  fi¬ 
nite  rate  equations.  To  illustrate  this  point.  Tables 
3,  4,  and  5  compare  the  relative  importance  of  effects 
in  the  electron  energy,  electron  density,  and  atomic 
hydrogen  density  equations,  respectively.  Each  table 
lists  the  ratio  of  every  term  in  the  equation  to  the 
dominant  term  at  each  grid  location,  averaged  over 
all  grid  points.  The  flow  is  divided  into  two  regions 
for  this  purpose,  the  arc  core  (a  >  0.01)  and  the 
outer  flow  (a  <  0.01), 

With  respect  to  the  electron  energy  equation.  Ta¬ 
ble  3  shows  that  ohmic  dissipation  is  an  impor¬ 
tant  source  of  energy  both  inside  and  outside  the 
arc.  This  leads  to  electron  temperatures  as  high  as 
20,  000° K  in  the  anode  attachment  zone  of  the  outer 
flow,  much  higher  than  the  1000  —  2000" /C  temper- 


Table  5:  Average  Ratio  of  Terms  to  Dominant  Term 
in  Atomic  Hydrogen  Density  Equation _ 


TERM 

ARC  I  OUTER  i 

Radial  Divergence 

0.545  i  0.664  ' 

Axial  Divergence 

0.487  I  0.910 

Radial  Diffusion 

0.265 1  0.001  ; 

Axial  Diffusion 

0.016  i  OTOOO  ! 

Net  Dissociation  (ff  -  ^2) 

0.131  I  OI9  ! 

Net  Dissociation  (e  —  H2) 

0.157  !  0.048  * 

Net  Ionization 

0.332  I  0T39  ' 

atures  which  would  be  calculated  by  a  model  with 
only  one  energy  equation.  Figure  11  illustrates  this 
result  by  comparing  radial  profiles  of  the  electron 
and  heavy  species  temperatures  at  an  axial  location 
0.25mm  downstream  of  the  constrictor  exit.  This  el¬ 
evated  electron  temperature  then  produces  enough 
electron  impact  dissociation  and  ionization  to  cre¬ 
ate  the  necessary  charge  carriers  for  electrical  con¬ 
duction  between  the  outer  arc  boundary  and  the  an¬ 
ode.  Tables  4  and  5  show  the  importance  of  these 
nonequilibrium  dissociation  and  ionisation  processes 
in  the  electron  and  atomic  hydrogen  continuity  equa¬ 
tions.  Radial  ambipolaz  diffusion  also  plays  a  role 
in  moving  ions  and  electrons  outward  from  the  arc 
into  the  surrounding  cooler  gas  Sow.  This  process  is 
evidenced  by  radial  profiles  of  the  ambipoiar  diffu¬ 
sion  and  net  ionisation  terms  in  the  electron  density 
equation,  shown  in  Figure  12  near  the  anode  wall  in 
the  current  attachment  region. 

Several  other  cases  have  been  run  to  evaluate 
the  model’s  effectiveness  in  predicting  arcjet  perfor¬ 
mance  over  a  range  of  parameters.  Figure  13  com¬ 
pares  simulation  predictions  of  specific  impulse  for 
the  TTl  thruster  at  three  different  applied  currents 
for  a  mass  flow  rate  of  O.lg/s.  Predicted  specific 
impulse  is  approximately  5-10%  higher  than  experi¬ 
mental  data,  but  should  come  into  closet  agreement 
when  the  results  of  the  aforementioned  anode  heat 
balance  model  are  incorporated  into  the  arcjet  flow 
model.  Table  6  compares  predicted  discharge  volt¬ 
ages  with  experimental  measurements  for  the  three 
cases  presented  in  Figure  13.  Voltage  predictions 
from  the  arcjet  simulation  consistently  fall  within  2% 
of  experimental  results,  demonstrating  the  code’s  ef¬ 
fectiveness  in  consistently  modeling  the  electric  arc 
characteristics  of  this  arcjet  thruster.  The  negative 
slope  of  the  V-I  relationship  is  also  captured. 

5  Conclusions 

A  detailed,  multifluid,  viscous  model  has  been  de¬ 
veloped  to  simulate  the  nonequilibrium  gasdynamic 
flow  in  an  electrothermal  arcjet.  The  model  is  im- 
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Table  6:  Comparison  of  Discharge  Voltages  for  the 
Cases  in  Figure  8 


Case 

Numerical 

Experiment 

I  =  60j4.  rh  =  O.ly/s 

117 

117 

I  =  lOOA.  m  =  O.lg/  s 

114 

_ : 

I  =  130A,  rh  =  O.lg/s 

112 

110 

plemented  on  a  nonuniform  mesh  fixed  to  experi¬ 
mental  thruster  dimensions,  and  the  equations  are 
solved  using  MacCormack’s  method  and  successive 
over-relaxation.  Numerical  results  are  achieved  for 
hydrogen  propellant,  and  calculated  thrust,  spe¬ 
cific  impube.  and  dbcharge  voltage  compare  well 
with  experimental  data  for  the  German  TTl  arcjet 
thruster.  The  internal  two-dimensional  structure  of 
the  flow  is  revealed,  particularly  with  iespect  to  arc 
development  and  anode  attachment,  and  the  two- 
temperature  nature  of  the  flow  b  evident.  In  par¬ 
ticular,  the  integration  of  a  separate  electron  energy 
equation  has  shown  that  stable  attachment  of  the  arc 
to  the  anode  occurs  by  increased  local  ohmic  heating 
coupled  with  nonequilibrium  dbsociation  and  ioniza¬ 
tion  in  the  flow  between  the  arc  core  and  the  anode 
wall. 

Additional  cases  are  being  run  at  thb  time  to 
obtain  more  data  points  for  comparbon  to  experi¬ 
ment  for  a  variety  of  operating  points  and  thruster 
configurations.  Abo,  an  anode  heat  balance  model 
has  been  developed  to  more  accurately  calculate  the 
anode  wall  temperature  for  thb  class  of  radiation- 
cooled  arcjet  thrusters.  Thb  model  b  currently  be¬ 
ing  integrated  with  the  interior  flow  model  in  order 
to  produce  more  fully  consbtent  solutions.  The  sim¬ 
ulation  may  then  be  used  to  identify  ways  to  improve 
overall  arcjet  thruster  performance. 

Acknowledgements 

Thb  research  was  supported  by  a  National  De¬ 
fense  Science  and  Engineering  Graduate  Fellowship 
adminbtered  by  the  U.S.  Army  Research  Office. 

References 

T]  R.  Spurrett  and  R.A.  Bond,  “Modelling  Ar¬ 
cjet  Thruster  Performance”,  lEPC  91-110, 
AIDAA/ALAA/DGLR/JSASS  22nd  Interna¬ 
tional  Electric  Propubion  Conference,  Viareg- 
gio,  Italy,  October  1991. 

2]  B.  Glocker,  H.O.  Schrade,  and  P  C.  Sleziona, 
“Numerical  Prediction  of  Arcjet  Performance”, 
AIAA  90-2612,  AIAA/DGLR/JSASS  21st  In¬ 
ternational  Electric  Propubion  Conference.  Or¬ 
lando,  July  1990. 


[3]  B.  Glocker  and  M.  Auweter-Kurtz.  "Numeri¬ 
cal  and  Experimental  Constrictor  Flow  Analy- 
sb  of  a  lOkW  Thermal  Arcjet”,  AIAA  92-3835, 
AlAA/SAE/ASME/ASEE  28th  Joint  Propul¬ 
sion  Conference,  Nashville.  July  1992. 

[4]  A.  Sakamoto,  Simplified  Modeling  and  Analysts 
of  Arcjet  Thrusters,  S.M.  Thesb,  Massachusetts 
Institute  of  Technology,  February  1993. 

[5]  M.  Andrenucci  et  al.,  “Development  of  a 
Computer  Programme  for  the  Analysb  of 
Arcjet  Nozzles",  lEPC  91-113,  AIDAA/- 
AlAA/DGLR/JSASS  22nd  International  Elec¬ 
tric  Propubion  Conference,  Viateggio,  Italy, 
October  1991. 

[6j  G.W.  Butler  and  D.Q.  King,  “Single  and 
Two  Fluid  Simulations  of  Arcjet  Performance”, 
AIAA  92-3104,  AlAA/SAE/ASME/ASEE  28th 
Joint  Propubion  Conference,  Nashville,  July 
1992. 

[7]  R.  Rhodes  and  D.  Keefer,  “Modeling  Arc¬ 
jet  Space  Thrusters”,  AIAA  91-1994,  AlAA/- 
SAE/ASME/ASEE  27th  Joint  Propulsion  Con¬ 
ference,  Sacramento,  June  1991. 

[8]  R.  Rhodes  and  D.  Keefer,  “Comparison 
of  Model  Calculations  With  Experimental 
Data  From  Hydrogen  Arcjets”,  IEPC-91-111, 
AIDAA/AUA/DGLR/JSASS  22nd  Intemar 
tional  Electric  Propulsion  Conference,  Viareg- 
gio,  Italy,  October  1991. 

[9]  R.J.  Biasca,  Chemical  Kinetics  of  Scramjet 
Propulsion,  S.M.  Thesb,  Massachusetts  Insti¬ 
tute  of  Technology,  September  1988. 

[10]  R.C.  Rogers  and  C.J.  Schexnayder  Jr.,  “Chem¬ 
ical  Kinetic  Analysb  of  Hydrogen-Au  Ignition 
and  Reaction  Times”,  NASA  Technical  Paper 
1856,  1981. 

[llj  R.K.  Janev,  W.D.  Langer,  K.  Evaiu,  and 
D.E.  Post,  Elementary  Processes  in  Hydrogen- 
Helium  Plasmas,  Springer- Verlag,  New  York, 
1987. 

[12]  M.  Mitchncr  and  C.  Kruger,  Partially  Ionized 
Gases,  John  Wiley  and  Sons,  New  York,  1973. 

[13]  E.J.  Sheppard,  Nonequilibrium  Ionization  in 
Electromagnetic  Accelerators,  Ph.D.  Thesb, 
Massachusetts  Institute  of  Technology,  1993. 

[14]  C.F.  Curtbs  and  J.O.  Hbschfelder,  “Transport 
Properties  of  Multicomponent  Gas  Mixtures”, 
Journal  of  Chemical  Physics,  Vol.l7,  June  1949, 
p. 550-5. 


9 


[15]  N  T.  Grier,  “Calculation  of  Transport  Proper¬ 
ties  of  Ionizing  Atomic  Hydrogen",  NASA  TN 
D-3186,  AprU  1966. 

[16]  V.A.  Belov,  “Viscosity  of  Partially  Ionized  Hy¬ 
drogen",  High  Temperature,  Vol.5,  1967,  p.31- 
6. 

[17]  J.T.  Vanderslice,  S.  Weissman,  E.A.  Mason, 
and  R.J.  Fallon,  "High- Temperature  Transport 
Properties  of  Dissociating  Hydrogen”,  Physics 
of  Fluids,  Vol.5,  No.2,  February  1962,  p.155-64. 

[18]  S.A.  Miller  and  M.  Martinc.;-Sanchez,  “Vis 
cous  and  Diffusive  Effects  in  Electrothermal 
and  MPD  Arcjet  Thrusters”,  lEPC  91-060, 
AIDAA/AIAA/DGLR/JSASS  22nd  Interna¬ 
tional  Electric  Propulsion  Conference,  Viareg- 
gio,  Italy,  October  1991. 

[19]  R.W.  MacCormack  and  B.S.  Baldwin,  “A  Nu¬ 
merical  Method  for  Solving  the  Navier-Stokes 
Equations  with  Application  to  Shock- Boundary 
Layer  Interactions”,  AIAA  75-1,  AIAA  13th 
Aerospace  Sciences  Meeting,  Pasadena,  Jan¬ 
uary  1975. 


[20]  P.  Kutler,  L.  Sakell,  and  G.  Aiello,  “On  the 
Shock-On-Shock  Interaction  Problem",  AIAA 
74-524,  AIAA  7th  Fluid  and  Plasma  Dynamics 
Conference,  Palo  Alto,  June  1974. 

[21]  B.  Glocker  and  M.  Auweter-Kurtz.  “Radia¬ 
tion  Cooled  Medium  Power  Arcjet  Experi¬ 
ments  and  Thermal  Analysis”,  AIAA  92-3834, 
AIAA/SAE/ASME/ ASEE  28th  Joint  Propul¬ 
sion  Conference,  Nashville,  July  1992. 

[22]  J.J.  Wang,  Electrodynamic  Interaction*  Be~ 
tween  Charged  Space  System*  and  the  Iono¬ 
spheric  Plasma  Environment,  Doctoral  The¬ 
sis,  Massachusetts  Institute  of  Technology,  June 
1991. 


Figure  1;  Electrothermal  Arcjet  Diagram 
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Figttie  2:  Physical  Grid  for  Atcjet  Simolations 


Figure  3:  Electric  Potential  Axial  Profile  for  the  Baseline  Case  Arcjet  Simulation 
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Figure  4:  Enclosed  Current  Contours  for  the  Baseline  Case  of  J  =  lOOA,  m 
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Figure  3;  Ionization  Fraction  Contours  for  the 
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Figure  12;  Radial  Profiles  of  Some  Terms  in  the  Electron  Density  Equation  in  the  Current  Attachment 
Region  Just  Downstream  of  the  Constrictor  Exit 
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Figure  13:  Predicted  Isp  Compared  to  Experimental  Data,  German  TTl  Radiation-Cooled  Hydrogen  Arcjet 
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Abstract 

A  numerical  code  for  the  plasma  properties  in 
the  acceleration  channel  of  a  Hall  thruster  was 
written  and  its  results  compared  with  experiments. 
The  model  was  constructed  to  help  understand  the 
acceleration  process  in  the  channel  and  to  relate  the 
axial  distributions  to  thruster  performance. 
Specific  to  this  analysis  are  the  definition  of  an  ion 
distribution  function,  careful  treaunent  of  the 
electron  energy  equation,  and  momentum  transfer 
between  neuuals  and  ions. 

Nomenclature 

(SI  units  unless  noted  otherwise) 

A  Channel  ctoss  sectional  area 
a  Ionization  fraction 
B  Magnetic  field  strength 
Cj  Random  velocity  of  species  s 
c.  Thermal  velocity  of  species  s 
D  Diffusion  coefficient 
5  Secondary  electron  emission  from  the  wall 
At  Time  step 

Az  Spatial  step  in  the  axial  direction 
E  Electric  field  strength 
e  Electric  charge  of  proton 
f ,  Ionization  energy 
17  Thrust  efficiency 
i7„  Acceleration  efficiency 
17,  Nonunifoimity  factor 
n.  Utilization  efficiency 
fs  Distribution  function  for  species  s 
F  Thrust 

F,  External  force  p.u.  volume  on  species  s 

r,  Spitzer  logarithm  for  species  s 

Y  Ion  production  cost 

fa  Anode  current 

hp  Specific  impulse 

Js  Elecuic  current  for  species  s 

js  Electric  current  density  for  species  s 
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k  Boltzmann  constant 
Kg  Heat  transfer  coefficient 
L  Channel  length 
m  Mass  flow  rate 
m,  Ion  mass  flow  rate 
rttg  Electron  mass 
m|  Ion  mass 

ns  Particle  density  for  species  r 
Ionization  rate 
ri^  Wall  loss  rate 

v„  Collision  freq.  of  species  s  with  species  r 
0  Electric  potential 
Anode  potential 

Qin  Ion-neutral  collision  cross  section 
qs  Electric  charge  of  species  s 
q.  Thermal  conduction  vector  for  species  r 
$2  Source  term  for  energy  equatk» 
a,  Ionization  cross  section 
Tno  Neutral  inlet  temperature 
T s  Temperanire  of  species  s 
Bohm  velocity 
V,  Velocity  for  species  s 
V,  Mean  mass  velocity  for  species  s 
W  Channel  width 

Number  of  electrons  in  kth  level 

(  )  Average  value 

Introduction 

Hall  thrusters  are  a  type  of  coaxial  plasma 
accelerators.  A  radial  magnetic  field  and  an  axial 
voltage  drop  between  the  anode  and  a  cathode  are 
externally  applied.  (See  Figure  1)  Neutral  gas 
(typically  Xenon  or  Argon)  flows  in  through  slits 
in  the  anode.  Ionization  takes  place  by  electron- 
neuuai  collisions  and  the  ions  are  then  accelerated 
axially  by  the  electric  potential.  Because  ionization 
takes  place  throughout  the  acceleration  channel, 
the  ions  accelerate  to  different  exit  velocities 
depending  upon  the  potential  difference  between  the 
location  of  their  birth  and  the  cathode.  A  cathode 
downstream  of  the  accelerator  releases  electrons  into 
the  ion  stream  neutralizing  the  flow.  A  fracdoo  of 
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these  electrons  travel  both  azimuthally  due  to  the 

ExB  Helds  and  drift  axially  toward  the  anode. 

The  numerical  code  developed  here  solves  a  set 
of  differential  equations  derived  from  taking 
moments  of  the  Boltzmann  equation.  Ohm's  Law 
and  an  ion  distribution  function.  For  the  code,  an 
ion  distribution  function  is  formulated  and 
integrated  to  generate  axial  ion  density  and  mean 
velocity  profiles.  Time  dependent  equations  for 
neutral  continuity  and  electron  temperature  are 
formulated.  Addidonal  equations  for  axial  current 
conservation  and  electric  Held  are  used  in  the  code. 

The  acceleration  channel  is  discretized  into  a 
fixed  grid  of  differential  .segments.  Axial 
distributions  are  initialized  onto  the  grid  at  the 
beginning  of  the  solution  routine.  A  convective- 
diffusive,  finite  difference  scheme  is  used  to 
integrate  the  set  of  equations  in  time  until  they 
converge  to  a  steady  state  solution.  From  the 
numerical  solution,  the  thrust  and  efficiency  of  the 
thruster  may  be  calculated. 

Assumptions  used  in  the  code  include  plasma 
quasi-neutrality  throughout  the  accelerator  and 
ambipolar  loss  of  electron-ion  pairs  to  the  walls  of 
the  channel.  Additionally,  in  this  one  dimensional 
study,  the  plasma  is  assumed  uniform  in  the  radial 
and  azimuthal  directions. 

For  the  cases  studied,  the  geometry  and 
operational  conditions  of  an  experimental  thruster 
are  reproduced  in  the  code.  The  results  produced 
from  the  numerical  analysis  are  then  compared  with 
the  experimental  results  from  research  performed  in 
Japan.  [1] 

This  paper  demonstrates  that  the  one- 
dimensional  numerical  code  can  be  a  useful  design 
and  plasma  analysis  tool.  The  method  is  versatile 
and  simple  enough  that  multiple  analysis  on 
various  thruster  parameters  (different  anode  currents, 
mass  flows,  etc.)  or  different  geometries  can  be 
examined  in  a  relatively  short  time. 
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The  model  for  the  thruster  includes  seven 
primary  equations  describing  the  unknown  variables 
-  neutral  and  electron  (or  ion)  density,  neutral,  ion 
and  electron  velocity,  electron  temperature  and 
elecuic  field.  From  these  values,  the  performance 
characteristics  of  the  thruster  may  be  calculated. 

Several  assumptions  are  used  to  facilitate  the 
computations  in  this  analysis.  The  first  is  (hat 
quasi  one  dimensional  flow  exists  resulting  in  no 
variations  in  the  radial  or  azimuthal  directions. 
Additionally,  it  is  assumed  that  only  one 
component  of  electric  and  magnetic  fields  exists  as 
illustrated  in  Figure  1.  Finally,  (he  flow  is 


assumed  to  be  quasi  neutral  in  the  acceleration 
channel. 


Figure  1 :  Hall  Thruster  Diagram 


Neutral  Continuity 

The  continuity  equation  for  neutrals  is  derived 
by  taking  the  first  moment  of  Boltzmann's 
equation.  The  resulting  expression  in  differential 
form  is 


dt  dz 


(1) 


The  terms  on  the  right  represent  the  rate  per  unit 
volume  at  which  ions  are  lost  to  the  walk  and  the 
ionization  rate  respectively. 

The  ion  loss  rate  to  the  wail  is  calculated  fmn 
the  rate  at  which  ions  enter  the  insulating  wall 
sheath: 


W 


S{z) 


where  the  Bohm  velocity  is  given  as 


(2) 


(3) 


In  Equation  2,  S(z)  represents  a  shape  factor  to 
account  for  imperfect  wall-plasma  contact  near  the 
injector  and  has  values  ranging  from  zero  near  the 
anode  to  one  about  halfway  down  the  channel  and 
beyond. 

The  ionization  rate  is  derived  assuming  a 
Maxwellian  electron  distribution,  a  nonelastic 
ionization  cross  section  according  to  Drawin's 
theory  [2],  and  integrating  over  all  electron 
energies.  The  resulting  expressions  used  are 
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/(©)=  (e  ^ug(u)du=  [e  ^  - — ln{l  25 B,u)du (5) 
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where  fi  is  the  ion  distribution  function  and  tbe 
Dirac  delta  function.  5,  is  employed  so  that  ions 
are  created  at  the  neutral  velocity.  In  Equation  1 1. 
n,  is  the  ionization  rate  and  h.  is  the  ion  loss  rate 
to  the  channel  walls  per  unit  volume.  This 
equation  however,  does  not  include  the  effects  of 
ion-neutral  collisions.  Adding  in  a  term  to 
approximate  this  transfer  results  in 

df,  eE  df,  .  .h., 

dt  dz  m,  n, 

-  w.)<5(Ui  -  u.)-(u,-u.)/,]  (12) 


a,  is  a  constant  obtained  by  gathering  all  other 
constants  into  one  and  has  the  value  for  Argon  of 
5  =  1.04x10*^^.  The  quantity  ()  is  evaluated  as 

()/\r=2.77xl0'^^  for  Argon,  f ,  is  the  ionization 
energy,  and  B2  are  constants  of  order  unity 
introduced  from  Drawin’s  cross-sectional  theory, 
and  u  represents  the  ratio  of  energy  to  the 
propellants  ionization  energy. 


The  equation  for  electron  momentum  is 
obtained  by  talcing  the  momentum  moment  of 
Boltzmann's  equations.  In  applying  the  result  to 
electrons,  the  terms  due  to  electron  inertia  and  the 
acceleration  of  heated  electrons  may  be  neglected. 
Additionally,  variations  in  the  radial  and  azimuthal 
directions  are  ignored  and  Bohm  diffusion  is  used  to 
obtain  an  effective  electron  collision  frequency. 
The  resulting  expression  after  these  simplifications 


nk^  +  kT.^  =  -en,E-l6eBn,v,  (9) 
dz  dz 


Rewriting  £  as  — ^  and  rearranging  Equation  9. 
dz 

an  expression  for  the  electric  field  is  obtained. 

^  =  16Bi;,+-^%-  +  -^  (10) 

dz  e  n,  dz  e  dz 


To  determine  the  ion  density  and  velocity,  a 
distribution  function  for  ions  is  derived.  In  the 
absence  of  collisions,  the  Boltzmann  Equation  is 
used  to  produce 


The  last  term  of  Equation  12  is  derived  in  such  a 
way  that  it  ensures  conservation  of  ions  and  gives 
the  correct  momentum  exchange  rate.  The  addition 
of  the  factor  of  1/2  preceding  the  last  term  is 
necessary  to  ensure  the  correct  momentum  transfer. 
Additioiiaily.  a  constant  value  for  the  ion-neutral 
collision  cross  section,  Qin.  is  used.  Rearranging 
slightly  produces 

dft  eE  dft 
dl  dz  m,  aX), 

[w,  -  vj 

-  f,  (13) 

2  J 

Equation  13  may  be  solved  using  the  method 
of  characteristics.  The  first  step  is  to  And  the 
characteristic  trajectories  and  then  to  calculate  the 
rate  of  change  of  fi  along  them.  The  characteristic 
system  is 

^  _  dUj  _ _ 

1  '  u.  ~  IT  ~  TTT  .  ( ti^  YT 

(n.  +  n.vJS(v,-vj-  ~+ v*.  f, 

'  \  n.  j 


Vi,  = -nMu,{^i  - '^n) 
v„  =  tJ,  -uj 
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Notice  that  the  first  equality  in  Equation  14  simply 
relates  position  and  velocity  and  the  second  and 
third  terms  produce  the  energy  equation  for  ions. 
However,  the  equality  of  interest  is  between  the  last 
two  terms: 


+  n,  )5(  v,-vj- 


(n  ) 

■ 

V 

I".  J 

/. 

I  If 

m, 

(17) 


This  equation  may  be  integrated  to  create  a 
distribution  function.  The  integration  is  over  the 
interval  from  =  c, ,  representing  ions  that  were 
just  created,  to  the  maximum  velocity  the  ions 
could  achieve  at  their  present  location.  This 
maximum  velocity  would  equal  the  ion  velocity  if 
they  were  created  at  the  anode  and  had  no  collisions; 


i».  =  -— (0, +  (18) 

\  /n, 


where  is  the  anode  potential  and  <p.  is  the 
potential  at  the  present  axial  location.  c„  is  the 


neutral  speed  of  sound  at  the  anode.  = 


so  that  The  distribution  function  is 

solved  for  each  panicular  location  z.  For  each  u, . 
there  is  a  different  Zo.  or  location  where  the  ion  was 
created,  such  that 


0{z,)=<P(z)->-^{v:  -  v;^)  (19) 


The  solution  to  Equation  1 7  is 


'n,(n,  +  n,vj 
-E 


e  (20) 


In  the  integral  appearing  in  Equation  20.  the 

quantities  £,  —  and  v,„  are  evaluated  at  the  time 

n. 

/'  and  location  z'  at  which  an  ion  with  velocity 
V,  at  (c.  t)  went  through  velocity  v' .  Thus,  u' 
acts  as  the  dummy  variable  of  integration  and  as  a 
parameter  specifying  the  integrand.  The  time  and 
location  (Zq.  to)  simply  the  ( /')  at  which 
this  particular  ion  was  created  and  at  which 
0,'  =  u, .  Also,  note  that  the  neutral  speed  tJ,  is 
to  be  evaluated  at  (;',  Das  well. 

The  ion  distribution  function,  integrated  over 
veliKity.  produces  the  ion  or  electron  density 


(21) 


The  mass  average  ion  velocity  can  be  found 
by  multiplying  the  right  side  of  Equation  20  by  u 
and  integraung 


n,  vjz) 


(22) 


Neutral  Velocity 

Due  to  ion-neutral  collisions, 
momentum  is  transferred  from  the  ions  to  the 
neutrals  resulting  in  an  increase  in  the  neutral 
velocity  and  a  decrease  in  the  ion  velocity.  For 
simplicity,  the  neutral  velocities  are  all  lumped 
together  into  their  mean  velocity  .  u, ,  so  that  a 
neutral  distribution  function  does  not  need  to  be 
created.  The  neutral  momentum  conservation 
equation  may  be  written  as 

"•Ki 

Mjn  is  the  momentum  change  due  to  collisions 
and,  assuming  the  same  mass  for  both  neutrals  and 
ions,  can  be  defuied  as 


0  ^ 

-  v,)f,]dv,  (24) 

Here  Qjn  (iie  collision  cross  section  and  for 
Argon  Q,>2=F4x10'^^  m^.  Equation  24  satisfies 
the  zeroth  moment  of  the  Boltzmann  equauon  and 
produces  the  correct  answer  for  the  first 
(momentum)  moment,  but  is  invalid  at  higher 
moments.  Simplifying  and  solving  EquaUoa  24 
yields 


^ n„Q„  j(  -  u,  y  f.dv,  (25) 

^  0 

and  evaluating  the  integral 

(26) 

The  term  within  <  >  represents  the  averaged  value 
of  the  quantity  inside.  Thus,  Equation  23  may  be 
written  as  _ 

vAz)  =  ^c;  +\,{z)  (27) 
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where  A,  is  given  as 

A,  =C>,  <28) 

j 

Electron  Energy  Enuation 

A  total  energy  equation  for  electrons  may  be 
obtained  by  taking  the  energy  moment  of 
Boltzmann's  equation  and  adding  it  to  the  dot 
product  of  the  electron  velocity  and  electron 
momentum.  The  resulting  expression 


The  term  Kg  is  the  heat  conduction  coefficient 
and  described  as  a  sum  of  two  limiting  models. 
The  first  is  due  to  collisions  and  is  the  “classical" 
description  in  the  presence  of  a  strong  magnetic 
field  [2], 


AT.,  =2.689x10-" (34) 

where 


£ 

di 


n.m. 


-~T 

2  J 


(36) 


n.m.  V. 


2m'  2 


(n,kT,d,)  =  S, 
(29) 


represents  the  total  electron  energy.  The  term  q, 
is  the  thermal  conduction  vector 


l=-KyT,  (30) 


The  alternative  heat  conduction  coefficient  is 
derived  according  to  Bohm  diffusion. 


'•  2jt  eB 


(37) 


Summing  these  two  thermal  conductivities 
produces  the  effective  heat  transfer  coefficienL 


The  right  side  source  term  includes  energy  lost  to 
collisions,  to  the  wall,  and  the  energy  gained  due  to 

the  electric  field  ( j  E).  This  equation  may  be 
simplified  by  introducing  the  variable  X 
representing  the  total  energy  density,  defined  as 


X  =  n.m 


f  3  A-  V 
2  m. 


:  3 


2  } 


'\n.kT. 


(31) 


So  that  Equation  29  may  be  rewritten  as 

[v,X^q,\^^  {n,kT,v,)^S,  (.32) 


Rearranging  and  expanding  equation  32.  its 
differential  form  is  obtained  as 


^.  =  ^.1+^.2  (38) 

In  Equation  38.  because  the  Hall  parameter.  P,  is 
large  K,  »  K,j  and  Bohm  diffusion  is  the  driving 
model  for  heat  conduction. 

The  source  terms  for  the  energy  equation 
include  energy  lost  to  the  wall,  energy  lost  in 
ionization  and  energy  imparted  to  the  plasma  by  the 
electric  field.  The  energy  input  from  the  electric 
field  is  given  as 


-  -  di 

j.  E  =  -en,v.E  =  en.v.— 


(39) 


Energy  lost  to  ionization  is 
(l+y)e£A 


(40) 


dt 


d[v.X) 

dz 


S,  -kT. 


dz 


(33) 


Here  the  term 


kT. 


dz 


has  been  moved  to  the 


right  hand  side  as  a  source  term  since  it  contains  no 
temperature  term. 


Where  the  term  (1  +  y)  represents  the  net  energy 

cost  for  each  ion  produced.  For  these  calculations, 
a  value  of  /  =  3  is  selected  based  on  a  more 
detailed  analysis  [3]. 

The  energy  lost  to  the  wall  consists  of  two 
elecuon  tluxes.  one  into  and  one  away  from  the 
wall,  as  shown  in  Figure  2.  The  first  term. 
Equation  41.  represents  the  energy  lost  by  the 
elecuon  flux  to  the  wall  and  the  second  is  the 
energy  flux  carried  by  the  secondary  electrons 
emitted  at  the  wall  back  into  the  plasma. 
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Summing  all  the  energy  sources  gives  the  loul 
energy  source  term  for  E^uauon  33 


'^plasma 


c  (,y  'T',  ^ 

•5; 


r,  =(i-5)r. 

Figure  2  Wall  fluxes 


=-r.(2e7-; 

(41) 

Tf,  -  r,e50. 

(42) 

(43) 

(44) 

and  S  is  the  secondary  emission  coefficient.  The 
total  energy  lost  to  the  wall  is  given  as 

Xr,,  =-er,[2rXi-5)0.]  (45) 


The  term  T'  represents  the  electron  temperature  in 
electron  volts,  and  0.  is  the  potential  drop  at  the 
wall,  given  by  T,  =(l-5)r,  and  T,  =n,Uj. 

which  is  valid  for  0„>T'  or  1-S>  2fr—. 

\  m, 

I'his  gives 


<0.  =T;in  (1-5)  — (46) 
4  V  Jtffi, 


Equation  41  and  42  are  represented  in  terms  of  the 
electron  flux  to  the  wall,  hut  by  using 

r,  =  (1  -  5)r,,  they  may  be  rewritten  m  terms  of 
the  tlux  of  ions  to  the  wall  or  the  ion  wall  loss 
defined  in  Equation  2: 

Z'  2T'  T 

Fj-  - -en^  — '—-*-0  1  (47) 

,  ^  1-d  J 


An  additional  relation  for  the  plasma 
properties  in  the  channel  is  one  for  current 
conservation  where 


-f  =  J,  +  j,  =  en,  V,  -  en,  u,  (49) 

A 

l(x  is  the  anode  current  and  ji  and  jg  ^  current 
densities.  The  equation  could  also  be  written  in  the 
form  /a  =  Jt  *  Je  lo  use  cur.ents.  The  ion  current 
J,  at  the  exit,  is  sometimes  referred  to  as  the  beam 
current  //>.  Figure  3  diagrams  the  currents  in  the 
thruster.  At  the  cathode,  a  part  lb  of  the  emitted 
electrons,  are  required  to  neutralize  the  current  of  Ib 
ions.  At  the  anode,  the  rest  of  the  cathode  electrons 
and  the  elecuons  produced  by  ionization,  lb .  equals 
the  anode  cuneni  Ja  ■ 


Anode 


I  Wc 


Cathode 


Figure  3  Hall  Thruster  Current  Diagiain 


The  thrust  produced  by  the  thruster  can  be 
calculated  knowing  the  exit  and  inlet  conditions. 
The  thrust  is  derived  as 

(50) 

The  specific  impulse  and  thrust  efficiency  are 
calculated  respectively  by 
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Additionai  boundary  conditions  at  tiie  end  of  tbe 
grid  are  that  the  electric  field  goes  to  zero,  the 
electron  temperature  equals  leV  and  the  slope  of  the 
electron  temperature  is  zero.  This  comprises  a  set 
of  boundary  conditions  for  six  first  order  and  one 
second  order  equation. 


Here  rn  is  the  mass  flow,  g  gravitational 
acceleration,  and  Iq  and  0^  the  acceleration  current 
and  potential. 

To  further  define  the  flow  efficiencies,  three 
internal  efficiencies  are  introduced  as  given  by 
Komurasaki  et.  al.  [1]. 


n.=  — 

m 


Accelerator 


Cathode 


Figure  4;  Thruster  Geometry 


-UlLL 

2e  [  m, 


17,  is  the  propellant  utilization  efficiency,  rj^  the 
primary  electron  utilization  and  17,  a  non 
uniformity  factor  representing  the  disuibution  in 
exit  velocities,  rj,  penalizes  the  acceleration  of 
ions  over  only  part  of  the  potential,  m,  is  the  ion 
mass  flow  rate  aixl  lb  the  beam  current. 

Subsequently,  by  their  definition,  the  thrust 
efficiency  may  be  expressed  as  the  product  of  the 
internal  efficiencies 

(56) 


The  numerical  model  consist  of  a  fixed  spatial 
grid  of  eighty  fixed  intervals.  Velocity  is  broken 
up  into  a  scaleable  grid  of  eighty  intervals  between 
tJ,  and  u,  .  The  neutral  continuity  equation 

(Equation  1)  and  tbe  electron  temperature  equation 
(Equation  33)  are  integrated  using  MacComiack’s 
method  [4],  a  second-order  accurate  predictor- 
corrector  method.  The  time  step  for  a  purely 
convective  equation  (Equation  I)  is  given  by  tbe 
Courant-Fieit^cbs-Lowy  condition  [4]  as 


The  time  step  for  a  conductive-diffusive  equation 
(Equation  33)  given  by  Anderson  is 


The  geometric  dimensions  of  the  thruster 
chosen  for  verification  of  the  model  [1]  are  shown 
in  Figure  4.  The  thruster  has  an  inner  radius  of  2 
cm.  an  outer  radius  of  2.4  cm.  and  accelerator 
length  of  8  mm.  The  flow  analysis  extends  out 
into  the  plume,  where  6  mm  past  the  accelerator 
channel  exit  plane,  the  cathode  is  located.  The 
accelerator  channel  has  a  constant  area  to  the 
physical  exit.  Beyond  the  exit,  ions  are  still  lost  to 
a  fictitious  wall  but  no  wall  recombination  takes 
place. 


The  inputs  to  the  code  include  the  mass  fiow, 
magnetic  field,  anode  current,  neutral  inlet 
temperature,  and  an  initial  or  pre-ionization  factor. 
From  these  parameters,  the  ion  and  neuual  deasities 
and  ion,  electron  and  neutral  velocities  at  the  antxle 
may  be  computed  and  used  as  boundary  conditions. 


\kn,v~- 


The  ratio  of  these  two  time  steps  A/^/A/^^N  can 
be  on  the  order  of  tens  of  thousands.  For  this 
reason,  the  election  energy  equation  is  integrated  N 
times  for  each  time  the  other  plasma  equations  are 
integrated  and  solved. 

The  flow  path  for  the  program  is  outlined  in 
Figure  5.  All  the  plasma  propertied  are  set  to  some 
initial  value  close  to  what  the  answer  is  believed  to 
be.  The  electric  field  is  then  frozen.  The  time 
steps  and  source  terms  are  evaluated  using  tbe 
values  of  the  previous  time  step.  Then  the  neutral 
conservation  and  momentum  equations  and  ion 
distribution  function  are  solved.  The  electron 
temperature  equation  is  integrated  N  times  and  tbe 
plasma  densities  are  checked  for  convergence.  If  tbe 
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densities  have  not  converged,  time  is  advanced  one 
step  and  the  equations  solved  again.  If  the  densities 
have  converged  then  the  electric  field  equation  is 
solved  and  checked  for  convergence.  If  it  has  not 
converged  the  electric  Held  is  frozen  at  its  new 
value  and  the  plasma  properties  solved  again.  If  the 
electric  field  has  converged  the  results  are  printed  to 
a  file. 


Figure  5:  Program  Flow  Path 


Results 

For  ail  the  runs  presented  here,  the  working 
gas  is  Argon  with  a  mass  flow  of  2.0Aeq 
(8.283x10'^  kg/s).  The  average  magnetic  field  is 
0.1  Tesla  with  a  minimum  value  of  0.093  T  and 
maximum  of  O.IOS.  The  loss  shape  faaor  used  in 
Equation  2  has  the  profile  shown  in  Figure  6  and  is 
the  same  for  all  the  numerical  runs.  The  secondary 
electron  emission  from  the  wall,  5 .  is  set  to  0.6 
and  a  value  of  3  is  used  for  / . 


Numerical/Expcrimental  Comparison 

For  this  run.  1.5 Amps  and  the  initial 
ionization  fraction  was  chosen  as  0.007.  The 
numerical  results  are  compared,  when  the  data  was 
given,  to  the  results  from  Komurasaki  where  the 
two  dimensional  profiles  have  been  averaged  over 
the  radial  direction  for  several  axial  locations. 


—  Numerical  Results  ^  Komurasaki  Data 


Figure  7:  Ion  Density 


Typically,  the  time  step  calculated  is  of  the 
order  At=»10’^  seconds,  except,  as  noted,  for  the 
electron  temperature  equation.  The  computational 
time  for  a  solution  is  on  the  order  of  10'^  to  lO'*^ 
seconds  and  the  time  for  an  ion  to  fiow  from  the 
anode  to  the  exit  is  of  the  order  of  10"^  seconds. 


Figures  7  through  1 1  show  the  profiles  of  the 
plasma  properties.  In  Figure  7,  the  ion  density 
rises  more  quickly  than  the  experimental  values. 
This  is  most  likely  due  to  incorrect  modeling  of  the 
shape  profile  near  the  anode,  in  that  there  are 
actually  more  losses  to  the  wall  than  modeled. 
However,  the  density  reaches  the  same  maximum 
density  that  the  data  does.  The  numerical  results 
peak  and  fall  off  earlier  than  the  data  which 
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indicates  that  eitber  die  ionization  rate  is  too  low  in 
this  region  or  the  wall  loss  is  too  large.  The  rapid 
falling  off  is  more  likely  due  to  the  ionization  rate, 
which  is  lower  than  expected  in  the  exit  region. 

xlO''23 


The  ionization  rate  (Figure  8)  reaches  its 
maximum  quickly  and  drops  away  smoothly  after 
that.  The  peak  is  reached  early  in  the  channel 
where,  even  though  the  ion  density  is  low. 
temperature  and  neutral  density  are  at  or  near  their 
maxima. 

Figure  9  shows  the  ion  and  neutral  velocity. 
The  neutral  velocity  increases  fairly  constantly  due 
to  the  ion-neutral  collisions.  The  ion  velocity 
shows  a  slope  discontinuity  at  the  channel  exit  (8 
mm  location).  This  is  due  to  continued  expansion 
and  acceleration  of  the  plasma  without  any  wall 
losses.  As  Figure  10  shows,  from  the  accelerator 
exit  to  the  end  of  the  computational  grid,  the 
potential  drops  less  than  half  of  the  anode  potential, 
while  the  ion  velocity  almost  doubles.  This  is  due 
to  the  majority  of  ions  actually  having  a  smaller 
velocity  than  their  mean.  These  ions  then  benefit 
more  in  the  remaining  potential  drop  tliat  those 
ions  with  velocities  greater  than  their  mean. 


Axial  Position  (mm) 

—  Ions  -  -  Neudals 
Figure  9;  Heavy  Particle  Velocity 

Figure  10  demonstrates  that  the  numerical 
code  accurately  calculates  the  anode  potential,  but 


that  the  potential  drops  away  slightly  quicker  than 
the  data.  This  indicates  that  over  the  acceler^ng 
channel,  the  electric  field  calculated  is  higher  than 
in  the  experiment,  while  it  is  smaller  outside  the 
channel. 


—  Numerical  Results  ^  Komurasaki  Data 
Figure  10;  Potential 


Figure  1 1:  Electron  Temperature 


The  electron  temperature  (Figure  11)  peaks 
near  the  anode  at  almost  fifteen  electron  volts.  The 
data  points  shown  in  this  case  represent  the  channel 
centerline  values  and  not  channel  averaged  values. 
The  experimental  temperature  at  the  anode  is 
slightly  lower  and  could  be  due  to  several  factors. 
The  first  is  that  the  computed  electric  field  is  too 
large  in  this  area  and  the  ionization  rate  and  wail 
losses  are  too  small.  Another  explanation  is  that 
the  value  selected  for  gamma  is  too  small  and  does 
not  accurately  account  for  secondary  ionization. 


Performance  Comparison 

The  performance  comparison  was  conducted  to 
check  the  numerical  code's  validity  at  predicting  the 
performance  characteristics  of  the  thruster  over  a 
wide  range  of  operating  powers.  For  this  analysis, 
the  average  magnetic  field  remained  at  0.1  Tesla, 
the  mass  flow  was  held  at  2.0  Aeq  while  the  anode 
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current  was  varied  between  0.75  and  2.5  Amps  to 
alter  the  operating  power  levels. 

Figure  12  shows  the  results  from  this  study. 
The  efficiency  was  found  to  vary  from  4.5%  at  an 
anode  current  of  0.75  Amps  to  almost  1 1%  at  2.5 
•Amps.  Over  this  same  range,  lae  anode  potential 
increases  from  127  Volts  to  225  Volts  and  the 
propellant  utilization  efficiency  rises  from  24%  to 
74%.  Both  the  acceleration  efficiency  and 
nonuniformity  factor  are  relatively  constant  over 
the  range  studied,  with  values  0  25  and  0.6 
respectively.  Their  insensitivity  to  the  accelerating 
current  agrees  with  results  and  conclusions  from 
Komurasaki's  analysis. 


^  Komurasaki  Data  —  Komurasaki  Data  Fit 
•  Numerical  Results 
Figure  12:  Performance  Results 


Conclusions 

The  purpose  of  this  study  was  to  create  a 
relatively  simple  yet  versatile  code  for  the  analysis 
of  a  Hall  thruster.  The  code  produced  demonstrates 
a  high  level  of  accuracy  in  both  predicting  the 
plasma  properties  in  the  acceleration  channel  and  in 
calculating  the  performance  characieri.stics  of  the 
thruster.  In  the  course  of  developing  the  code,  a 
model  was  constructed  of  the  ion  distribution 
function  accounting  for  momentum  transfer 
between  heavy  species  and  for  the  effect  of  wail 
losses.  Further  studies  will  test  the  code  over  a 
wider  range  of  operating  conditions  and  model  other 
thruster  geometries  to  test  its  full  viability. 
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Th«oretical  and  computational  predictions  of 
MPD  thruster  efficiency  are  substantially  higher 
than  experimental  measurements  for  comparable  ge¬ 
ometries.  One  reason  for  this  discrepancy  is  the  large 
anode  voltage  drops  found  in  experiments  but  not  in 
current  models.  This  paper  shows  that  these  voltage 
drops  could  be  occuring  in  the  quasi-neutral  near  an¬ 
ode  region  of  the  plasma.  These  near  anode  drops 
arue  because  of  starvation  of  the  plasma.  A  one 
dimensional  analytical  model  is  used  to  show  that 
these  voltage  drops  may,  in  principle,  exist.  A  com¬ 
putational  axisymmetric  thruster  simulation  is  used 
to  show  that  the  voltage  drops  which  arise  by  this 
mechanism  show  good  agreement  with  experimental 
data. 


B  .Magnetic  field  strength. 

C  Random  velocity, 

e  Charge  of  a  proton. 

E  Electric  field  strength. 

E,  Ionisation  energy. 

J  Current  density. 

k  Boltiman's  constant, 

m  Particle  mass. 

n  .Number  density. 

n.  Electron  production  rate. 

P  Pressure, 

f  Time. 

Q  Collision  cross  section. 

T  Temperature. 

V  Velocity 

I  Ionisation  fraction. 

1/  Collision  frequency. 

^  Mass  density. 

j  Electrical  conductivity. 
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<'upyriitht  1^1993  by  the  American  Institute  of  Aeronautics 
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0  Heat  conduction  coefficient. 

V  Viscosity  coefficient. 

Note:  The  subscripts  e.o.i,  and  g  refer  to  electrons, 
neutrals,  ions,  and  heavy  species  respectively.  When 
used  as  subscripts  or  otherwise,  r,  B,  and  t  refer  to 
the  coordinate  directions. 


Magnetoplasmadynamic  (MPD)  thrusters  have 
been  considered  as  a  means  of  propulsion  for  space¬ 
craft  on  long  missions  for  the  past  thirty  years.  MPD 
thrusters  produce  thrust  in  the  range  of  1-100  New¬ 
tons  at  specific  impulses  ranging  from  1000-5000  sec¬ 
onds.  They  are  expected  to  be  most  efficient  at 
power  levels  on  the  order  of  MegaWatts.  However, 
the  experimentally  measured  efficiency  of  these  de¬ 
vices  is  too  low  for  them  to  be  attractive  for  space 
use.  Experiments  show  efficiencies  under  30%  in  al¬ 
most  all  cases  [17,  10].  These  experimental  measure¬ 
ments  are  much  lower  than  efficiencies  predicted  by 
theoretical  or  computational  means  >'12,  14|.  As  the 
devices  are  run  at  fixed  currents,  the  efficiency  is  a 
function  only  of  the  Jet  power  and  the  voltage  be¬ 
tween  the  electrodes.  Computational  and  analyti¬ 
cal  predictions  of  jet  power  are  usually  quite  good 
(11,  14).  Therefore,  the  discrepancy  in  efficiency 
must  arise  due  to  low  voltage  predictions  by  the  var¬ 
ious  models.  Recent  experiments  have  shown  that  a 
large  component  of  the  voltage  drop  in  experimental 
devices  takes  place  very  near  the  anode  and  results 
mainly  in  heating  of  the  anode  [5,  7|. 

This  paper  attempts  to  show  that  these  voltage 
drops  may  be  a  result  of  starvation  of  the  near  anode 
plasma.  Section  2  details  the  mechanism  by  which 
these  drops  arise.  Section  3  describes  an  analytical 
model  which  gives  some  insight  into  the  physics  of 
the  voltage  drops.  Section  4  describes  a  computa¬ 
tional  model  of  MPD  thrusters.  Section  5  presents 
results  from  that  simulation  showing  good  agree¬ 
ment  between  numerically  predicted  and  experimen- 


tally  measuted  voltage  drops. 


2  Theory 


If  there  were  such  a  thing  as  an  ideal  thruster  with¬ 
out  the  Hall  effect,  there  would  be  no  axial  current 
and  the  current  would  run  straight  from  electrode 
to  electrode.  However,  due  to  the  Hall  effect  there 
is  axial  current,  the  magnitude  of  which  is  approxi¬ 
mately  given  by 


J.  »- 


<rB$ 

€n« 


=  I3J. 


(1) 


predicts  the  resulting  voltage  drop.  Bakhst  assumes 
that  the  axial  electric  field  is  zero  with  the  axial 
electron  pressure  gradient  and  ion  velocity  neglected. 
Assuming  that  the  ion  and  electron  temperatures  are 
constant  radially  and  combining  equations  1  and  2 
yields 

dr  en,k(T,+Tty  ’ 

If  Bt  and  Jr  are  constant  radially  then 


"•(»•) -t/  gk(Tr+T,) 


(5) 


where  0  =  is  the  Hall  parameter  and  the  mag¬ 
netic  field  is  into  the  plane  of  the  thruster  and  so 
has  a  negative  sign.  Since  the  Hall  parameter  is  in¬ 
versely  proportional  to  the  electron  number  density, 
0  and  Jt  increase  as  n«  decreases.  Increased  axial 
current  leads  to  higher  radial  Lorentz  force,  which 
must  be  approximately  balanced  by  a  radial  pressure 
gradient; 


dn,k(Tt  +  T,)  ,  „ 

- — - -  -J,Bt 


(2) 


From  equation  3,  the  contribution  to  the  radial  elec¬ 
tric  field  from  the  Hall  and  electron  pressure  terms 
is  then 


Er,B*U 


<rB}Jr  T, 
(en.)»T.-^T, 


(6) 


The  potential  drop  across  the  near  anode  layer  is 
then  given  by 


AV(anode) 


Erdt 


Since  J,  and  B$  are  both  negative  and  r,  and  Tf  are 
expected  to  be  relatively  constant,  the  radial  gra¬ 
dient  of  n,  must  be  negative.  Therefore,  n.  must 
decrease  dramatically  in  the  radial  direction  if  J,  is 
large,  as 

5n,  J,  Bt 

dr 

So,  low  electron  number  density  and  high  axial  cur¬ 
rent  can  build  on  each  other  and  lead  to  anode  star¬ 
vation. 

Once  the  anode  becomes  starved,  large  radial  elec¬ 
tric  fields  will  develop,  as  seen  from  Ohm’s  Law.  The 
Hall  component  of  the  radial  field,  which  is  dominant 
in  the  starved  region,  is  given  by 

fj,  «  (3) 

en,  or 

Since  n,  is  small  and  Jt  is  large,  the  radial  elec¬ 
tric  field  can  be  quite  large.  Starvation  has  been 
observed  experimentally  [9].  However,  the  width  of 
the  starved  region,  and  therefore  the  extent  of  the 
voltage  drop  across  it,  has  not  been  experimentally 
determined.  If  this  field  is  large  enough  and  occurs 
over  a  wide  enough  region,  the  voltage  drop  across 
the  starved  region  can  be  significant. 

3  Analytical  Model 


_  _ fc(T,  +  r,)en*p _ 

2e  {r^uT  -  r,)2(rBlJr  +  *(T,  +  T,)enJo ' 

(7) 

Bakhst  assumes  that  the  anode  current  density  will 
be  due  to  the  random  electron  thermal  flux  to  the 
wall,  so  that 

4/, 

Bakhst’s  model  can  be  used  to  predict  the  an¬ 
ode  voltage  drop  between  the  anode  and  the  last 
inner  point  of  the  simulation  described  in  Section 
5.  At  a  point  midway  down  the  electrode,  T*  = 
23,000K,T,  =  Z2WK,Bi  =  0.0726r,<r  =  3580|f, 
and  Jr  =  -2.08  x  10“*^.  So,  C  =  591,000m/s 
and,  to  get  the  right  thermal  flux,  n.o  =  8>8  x  10‘*. 
The  last  interior  point  is  0.1  mm  away  from  the  an¬ 
ode  at  r.  =  0.072m.  The  Bakhst  model  predicts  a 
potential  drop  over  the  last  cell  of  -1.09  Volts.  The 
distribution  of  the  electron  number  density,  the  ra¬ 
dial  electric  field,  and  the  potential  drop  are  shown 
in  Figures  1-3. 

However,  this  Bakhst  formulation  neglects  a  num¬ 
ber  of  effects,  particularly  in  the  transverse  momen¬ 
tum  equation.  A  somewhat  more  complete  model 
has  been  developed.  This  model  also  assumes  that 
the  axial  electric  field  is  zero  throughout  the  near 
anode  region, 

1  1  dBi  1  d 

— —  -  —  Vi9 B$  4*  ■  "T”  ^  “  4*  “ 

(8) 


Probably  the  first  model  of  anode  starvation  is 
that  of  Bakhst.  Bakhst’s  [2]  focus  was  on  deter¬ 
mining  when  starvation  occurs,  but  his  model  also 


where  electron-neutral  collisions  have  been  ne¬ 
glected.  The  transverse  momentum  equation  for  the 
charged  species  is  modified  by  including  the  contri¬ 
bution  due  to  ion-ion  viscosity,  5,>  so  that 


1  B$  drBj 

k{Tf+T,)  nor  dr 


Even  though  the  ion-neutral  viscosity  is  2-5  times 
larger  than  the  ion-ion  viscosity,  due  to  the  low  ion¬ 
ization  fraction,  the  radial  gradients  of  the  ion  veloc¬ 
ity  are  100-1000  times  larger  than  those  of  the  neu¬ 
tral  velocity.  Therefore,  the  viscosity  source  term, 
given  in  complete  form  in  Section  4  will  be  assumed 
to  be  given  by 

Over  the  small  transverse  region  being  modeled, 
the  change  in  the  axial  flow  must  be  relatively  small. 
Therefore, 

(10) 

Taking  the  appropriate  derivatives, 

d^Vir  _  ,  2Vir,dpi^ 

dr*  Pi  dr*  p?  dr 

The  radial  momentum  equation  can  then  be  writ- 


^  \b.j.  +  *(r,  +  r.)^l 

dr*  n,  dr  ’  4u<iV5,  [  '  '  '  dr  \ 

(11) 

where 

Bt  dB$  1  d  1 

=0-  VirBt+  h  • 

eritpo  dz  en,  dz  J 


with  the  Bakhst  results  in  Figures  1  -  3.  The  solu¬ 
tions  were  obtained  by  shooting  with  a  variable  step 
size  Runge-Kutta  algorithm.  The  total  Hall  poten¬ 
tial  drop  across  the  near  anode  region  predicted  by 
this  model  is  -9.68  Volts.  The  difference  between 
the  two  models  is  due  to  the  lower  electron  number 
density  over  the  layer  predicted  by  the  modified  the¬ 
ory.  The  modified  theory  predicts  lower  n,  because 
the  total  pressure  (fluid  and  magnetic)  away  from 
the  anode  is  less  than  that  at  the  anode  due  to  vis¬ 
cous  dissipation.  Lower  total  pressure  for  constant 
temperatures  implies  lower  electron  number  density. 
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A  three  fluid  simulation  of  MPD  thrusters  has 
been  developed  and  used  to  show  that  anode  star¬ 
vation  can  lead  to  voltage  drops  of  the  same  magni¬ 
tude  seen  experimentally.  The  simulation  assumes 
that  the  plasma  is  quasi- neutral.  The  geometry  is 
assumed  to  be  cylindrical,  with  no  variation  in  the 
azimuthal  direction.  The  magnetic  field  is  assumed 
to  be  confined  to  the  azimuthal  direction.  The  elec¬ 
tric  field,  the  current,  and  the  velocity  all  have  com¬ 
ponents  in  both  in  plane  directions.  The  plasma 
is  treated  as  being  out  of  equilibrium,  in  that  the 
electron  and  heavy  species  temperatures  are  treated 
separately  and  coupled  only  by  elastic  collisional  en¬ 
ergy  transfer,  the  ionization  Raction  is  not  deter¬ 
mined  by  Saha  equilibrium,  but  by  balancing  ion 
flow  with  ionization  and  recombination  collisions, 
and  the  ion  and  neutral  velocities  are  coupled  only 
by  collisional  drag.  The  goveriting  equations  and  the 
relevant  source  terms  are  described  below. 

4.1  Governing  Equations 


Computational  Model 


The  magnetic  field  is  given  by 


drB, 

dr 


=  PoJsr 


(12) 


and  the  Hall  component  of  the  radial  electric  field  is 
given  by 


= - BfJt  -f  kT, 


In  order  to  match  the  Bakhst  theory,  the  magnetic 
field  and  the  electron  number  density  at  the  anode 
were  chosen  as  two  of  the  boundary  conditions.  The 
third  boundary  condition  was  that  the  radial  deriva¬ 
tive  of  n*  was  zero  at  the  outer  edge  of  the  layer.  For 
the  same  T,,Tg  and  <t  as  used  for  the  Bakhst  model 
and  with  additional  conditions  that  vu  —  7.3  x  10~* 
and  Vir(d)  =  Vj  at  the  wall  and  assuming  ^  =  0, 
the  solutions  for  this  modified  theory  are  given  along 


Incorporating  the  above  assumptions  yields  a  model 
consisting  of  nine  partial  difFerential  equations. 
There  are  eight  differential  equations  for  the  fluid 
variables  and  one  for  the  magnetic  field. 

dn,  d{n,Vir)  ,  d[ntVi,)  n.K, 

IT  + 

(U) 

and 

dnn  d(nnV„)  d(n„v;,)  n„V„, 

(15) 

where  hn  represents  the  number  of  recombination 
events  and  hi  represents  the  number  of  ionization 
events.  As  mentioned  above,  the  plasma  is  assumed 
to  be  quasi-neutral,  so  that  n,  is  used  in  place  of  rii. 
It  is  useful  to  define  the  global  density, 

p  =  m,Ti,  -f  min,  +  %  m„{n,  +  n,). 


It  is  convenient  to  define  the  glob&l  velocity,  the  cur¬ 
rent,  and  the  neutral  slip  velocity,  where,  respec¬ 
tively, 

^  ^  m.n.V.-^^^n.V^■fm,nnVn  ^ 

P 


Un 

J  =  -en.(V.-V<) 
u  =  Vi  -  V,. 

The  state  equation  for  each  species  is  assumed  to  be 
given  by 

P.  =  n,kT..  (16) 

Therefore,  the  neutral  momentum  equations  are 

a(PnVn,)  d{fiM+Pn)  d(pnV^,Vn.) 

at  dr  dz 

=  Snn.  +  Sni.  +  K^i(Vu  -  K.,)  +  -  Vnr) 

o»V* 

-  n/winV;,  +  njtntnVi, - —  (17) 


,  S(PnVn.Vns)  .  OjPnV^,  +  P^) 

at  8r  dz 

=  SnnM  +  SnU  +  Kni{Vu  -  V,,)  -H  lf„.(V.,  -  V„) 

-  n/m„V;.  tiRtn^Vu  -  (18) 

where  Sntr  and  5nt<  represent  the  viscous  terms  and 

K,t  -  n,m,tv,t  =  Ku 

m,V  =  -(-  m,"' 


B(PiVu)  ,  Q{PiVM  ,  d(piV,\  Pi  +  P.)  ,  „ 

- Tr - ^ - fz - 

+5ii.  +  5i„  -t-  Kin(y^,  -  Vi,)  +  K,n(Vn,  -  V.,) 

-I-  njmiVn,  -  njrmiVi,  -  (22) 

If  Ti  =  Tn  then  the  overall  heavy  species  energy 
equation  is 

B(\P,)  ,  g(|P,vO  ,  g(fP,Vs)  ,  ..avr  ,  av,' 

at  dr  dz  dr  dz' 

=  ii(/i:i.  +  /i:„)(r.-r,)-(-*-i-^-^«,-5^ 

TTln  2  r 

-h  i  [K^i  +  Kin  +  mi(nr  +  nj,)]  {U*  +  Uf).  (23) 

where  T,  =  Ti  =  T,  and  P,  =  ;^pT,. 

Using  the  electron  continuity  equation  the  energy 
equation  can  be  rewritten  in  non-conservative  form 


2«.*l  ^  +  •*  dg  dr  ^  dz 

=  (jc.i  -b  K,n)  f^(r,  -  T.)]  +  «.  +  ~ 

.”*a  J  ^ 

-ni:.n(£r*-2— )-n.(jPi-|*r.)-^.  (24) 

€flf  m  r 

From  Maxwell's  equations, 

BB$  dEr  dEt 

■©r  =  "dT 

The  scalar  components  of  the  entrant  axe  given  by 


The  electron  momentum  equations  are  used  to  de¬ 
rive  the  Ohm’s  laws,  given  by 

E,  =  Vi,B,  -  — (J,P,  -(-  ^)  -^  (19) 

en,  dr  <r 

1 

=  -l^i,5,  -  — (^  -  JrB,)  +  ^  (20) 

en,  vz  <r 

These  equations  can  be  used  to  replace  the  electric 
field  term  in  the  ion  momentum  equation.  Doing 
this  yields  the  ambipolax  momentum  equations 

d(PiVi,) ,  d(piVi;  -f-  Pi  p.) ,  d(piVi, Vi,) 

- - - + — _ —  =  _j,s, 

-i-5ii,  +  5in,  lifin(V„,  -  Vir)  +  A’.n(V„,  -  V„) 

+  himiVnr  -  hftmiVir  -  (21) 


j  =  =  +  A  (27) 

*  dr  fto  dr  Mcr 

4.2  Source  Terms 

Because  of  the  substantial  slip  between  neutrals  and 
ions,  the  viscous  source  tenns  are  more  comidex  than 
those  in  the  one  fluid  Navier-Stokes  equations.  A 
derivation  and  expression  for  these  terms  was  found 
in  a  recent  paper  by  Fernandes  and  Fexnandea[3]. 
According  to  their  derivation  the  non-isotropic  part 
of  the  pressure  tensor  is  given  by 

n,  =  2v,i V  o  Vi  -I-  2v,»V  o  V„.  (28) 

where  the  operator  V  o  V  is  the  symmetric  traedeu 
gradient  of  V  given  by  [4] 

VoV,  =^, -i(V.V,)l  =  T.  (29) 


4 


Exp&nding  the  various  derivatives  yields  the  viscous 
source  terms  [1], 


c  -  r4  5*v.,  a*v;- 

“  '"*•  3  dr^  ^  dz^  3  drdz 


4dVtr.V,t  ^V4t  .  ^Vtr  Sv,t 
'^3  dr  r  ^  dr  dz  dz 

dVt,  dv,t  2  dVu  dv,t  ^V;,  5^ 

dr  dz  3  dz  dr  Z  r  dr  r 


*“  ’^'‘[3  dz*  dr*  '^3  drdz 

BV%^  .  1  t/jj  Bvgt ,  2  BVff  Bvgg 

"’’“dT^a”  ~  3  dr  dz 

2  Vtr  dVft  ,  dVtx  ,  I  t^<t  ^  ,  4  dVft  dVtt  /«.\ 

3  r  dz  dr  ^  dr  r  3  dz  dz 

The  pressure  tensor  also  appears  in  the  energy 
equations.  The  non-diagonal  terms  of  the  pressure 
tensor  are  buried  in  .he  viscous  dissipaton  term  and 
are  given  by  [1] 

♦  =  ♦<  +  ♦,  =  [(v,<  -I-  v.„)  (|((^)* 

^  r  ’  dz  ’  r  dr 

dv„dv.,  V„dv.,  dv,r  sv;.  ,\l 

dr  dz  r  dz  dz  dr  '  )\ 

(32) 

The  various  viscosity  coefficients  ate  given  by[3] 


J/2  .  .  ^  ,  <*(1  -a)v. 


Vii  =  +<*»)  + 


]/?<■  (33) 


v„n  =  [(1  -  «)^(|  +a,)  +  (34) 

V  1/f 


The  collision  integrals  O  are  given  by  [4] 


drdz 

J  'f'l _ _ 1 _ I  : _ _ I  .^11!.^ _  -11 


Win  =  v„i  =  a(l  -  a)(- -  Oi)/q<  (35) 


where 


,  (1  -  ay 


^  IX  — 

gi-(3+ai)[-  +  — 


vr. 


The  neutral- neutral  and  ion-neutral  collisions  will 
be  modeled  as  hard  sphere  collisions.  Using  exper¬ 
imental  values  from  Lieberman  and  Velikovich  [13], 
the  relevant  cross  sections  are  given  by 

=  ?(1.7  X  lO-i'Tn^)  (37) 


(2)  _  ^  ^ 


=  -O. 


=  -(1.4  X  10-^®).  {3«) 
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All  units  in  the  thesis  are  MKS  unless  otherwise 
noted.  For  the  ion-ion  collisions  the  integrals  can 
be  approximated  using  the  Coulomb  interaction  po¬ 
tential  cut  off  at  the  Debye  length.  From  Fersiger 
and  Kaper  again. 


where 


r,  =  1.24  X  lO^W^. 


The  heat  conduction  terms  k,  are  given  by  k,  = 
—V  •  Hf.  The  heat  flux  vector,  H  is 

H,  =  -e,VT..  (41) 

where  @  represents  the  thermal  conductivity.  So, 
the  contribution  from  the  heat  conduction  terms  is 

^  fi  dT,  ,  d*T,  .  d*T,l .  dT.  ae,  ,  ar.  ae, 

[r  dr  dr*  dz*  J'*'  dr  dv  dz  dz  ‘ 

(42) 

From  Mitcbner  and  Kruger  [15],  the  thermal  con¬ 
ductivity  for  the  electrons  is  given  by 

r\  3.4  k*n,T,  ^  , 

»•  -  <«' 

For  the  heavy  species  the  thermal  conductivity  is 
given  similarly  by 

e  ^  I  "*  ] 

'  TtliCf  [rinQnn  +  ritQin  ««<?«  + »tn<?tnj 

(44) 


where 


e^lnFi 

32xe*k*Jf 


and 


and 


Sheppard  [16]  has  developed  a  fit  to  a  multilevel 
model  which  compares  well  with  experiment.  His 
recombination  coefficient,  R,  for  Argon  is  given  by 
[16] 

R  =  8,25  X  io-«e*  lift-®  »»>’  (46) 

The  recombination  rate  is  then  given  by 
nji  =  Rn,. 

By  detailed  balancing,  the  ionisation  rate  is  then 
given  by 

hi  =  A5n«n„ 

where 

S  =  2.9  X  10»T.*  exp 

5  Computational  Results 

The  experimental  device  on  which  this  research 
concentrated  was  the  Constant  Area  Channel  (CAC) 
studied  by  Heimetdinger,  Kilfoyle,  and  Martines- 
Sanches  [7,  8,  6,  9].  However,  the  bulk  of  the  ex¬ 
perimental  work  reported  by  Heimeringer  et.  al.  in¬ 
volved  the  Fully  Flared  Cathode  (FFC),  which  dif¬ 
fered  from  the  CAC  in  the  shape  of  the  cathode, 
but  was  the  same  device  in  other  respects.  Results 
from  the  two  thrusters  will  be  treated  somewhat  in¬ 
terchangeably  in  the  following  discussion.  Where 
necessary,  a  distinction  will  be  drawn  between  the 
two  devices. 

The  physical  thruster  was  an  axisymmetric  device 
with  the  cathode  as  the  inner  electrode.  Both  the 
cathode  and  anode  were  0.09  m  long,  with  a  constant 
outer  radius  for  the  cathode  of  the  CAC  of  0.053  m 
and  a  constant  inner  radius  for  the  anode  of 0.072  m. 
The  FFC  had  a  cathode  which  varied  from  0.042  m 
outer  radius  at  the  inlet  to  0.053  m  at  the  throat  to 
0.033  m  at  the  thruster  exit.  Currents  ranged  from 
20  kA  to  60  kA  with  a  mass  flow  of  4  x  lO-®^  . 
Onset  appeared  to  occur  at  approximately  60  kA. 

The  numerical  thruster  is  depicted  in  Figure  4.  It 
consists  of  two  concentric  cylinders,  of  which  the  first 
0.109  m  of  each  is  conducting,  followed  by  a  short 
insulating  section  of  0.031  m.  The  interelectrode  gap 
is  0.02  m,  with  a  cathode  inner  radius  of  0.052  m. 
The  plume  is  not  included  in  the  simulation.  The 
mass  flow  of  4  xl0~®  kg/s  is  assumed  to  be  injected 
through  the  whole  backplate,  with  the  mass  flow  per 
unit  area  constant  at  all  radial  locations.  A  number 
of  different  current  leveb,  ranging  from  23.4  kA  up 
to  39.0  kA,  were  simulated.  These  currents  are  all 
well  below  the  onset  condition  but  span  the  region 
where  large  anode  voltage  drops  develop. 


As  discussed  in  Section  1,  simple  theories  and 
numerical  modeb  of  MPD  thrusters  have  always 
significantly  underpredicted  the  total  voltage  of 
thrusters  as  observed  in  experiments.  The  data  from 
Heimerdinger,  as  well  as  from  other  experiments  dis¬ 
cussed  previously,  indicate  that  the  total  voltage  b 
dbtributed  between  a  cathode  fall  voltage,  a  bulk 
plasma  voltage,  and  an  anode  fall  voltage.  None  of 
the  exbting  theoretical  or  numerical  results  show  the 
anode  and  cathode  fall  voltages,  or  adequately  ex¬ 
plain  the  cause  of  the  anode  falb.  Thu  b  why  the 
theoretical  and  numerical  results  predict  voltages  so 
much  lower  than  those  observed  experimentally. 

The  model  and  simulation  used  in  thu  research 
however  does  seem  to  reproduce  the  anode  voltage 
drop  behaviour  seen  by  Heimetdinger.  The  anode 
voltage  drops  seen  in  the  simulation  occur  in  the 
quasi-neutral  bulk  plasma  very  neat  the  anode.  The 
voltage  drops  appear  even  though  the  model  uses 
fluid  equations  and  does  not  include  the  non-neutral 
anode  sheath.  These  voltage  drops  occut  because  of 
the  basic  mechanbm  outlined  in  Section  3,  as  will 
be  shown  by  the  data  presented  herein. 

Figure  5  shows  the  cturent  lines  for  the  baseline 
case,  with  the  concentration  of  lines  a  measure  of 
the  current  density.  The  plot  u  drawn  so  that  the 
axial  and  radial  dimensions  are  roughly  in  the  cor¬ 
rect  proportions.  The  numbers  on  the  plot  indicate 
the  approximate  percentage  of  current  enclosed  by 
the  nearest  contour.  As  can  be  seen,  the  current  b 
highly  skewed  near  the  anode,  turning  almost  paral¬ 
lel  to  the  electrode.  TUs  u  because  the  axial  com¬ 
ponent  of  the  current  b  substantially  larger  than  the 
radial  component.  The  axial  component  b  so  large 
because  the  Hall  parameter  near  the  anode  b  quite 
large,  reaching  values  as  high  as  iOO. 

The  Hall  parameter  b  so  high  because  of  the  low 
dectron  number  density.  The  electron  number  den¬ 
sity  along  five  radial  cuts  b  shown  in  Figure  6.  The 
cuts  are  at  the  axial  locations  shown  in  the  figure 
key.  The  electron  number  density  decays  from  a 
maximum  of  2  x  10®®m~®  near  the  cathode  to  a 
minimum  of  1  x  at  the  anode.  The  elec¬ 

tron  number  density  b  dropping  because  both  the 
ionbation  fraction  and  the  total  mass  density  are 
dropping  near  the  anode.  The  ionbation  fraction 
b  low  near  both  the  cathode  and  the  anode  due 
to  recombination  of  the  ions  at  the  wall.  The  low 
ionbation  baction  near  the  cathode,  and  the  slight 
(relatively)  drop  in  electron  density  there  leads  to  a 
small  cathode  voltage  drop. 

The  electron  number  density  neat  the  anode  b  so 
low  because  the  plasma  has  been  turned  away  from 
the  anode  by  the  radial  Lorents  force.  Thb  turn  b 
shown  in  Figure  8  where  there  b  a  region  of  large 
negative  radial  velocity  near  the  backplate.  Num- 


bets  on  the  plot  show  values  at  the  local  minima 
('582, -177, -61, *346  m/sec)  and  the  maximum  value 
(754  m/sec).  It  is  also  illustrated  by  the  stream¬ 
lines,  shown  as  the  dotted  lines  in  Figure  10.  This 
turn  leads  to  depletion  of  the  plasma  near  the  an¬ 
ode,  as  shown  by  the  isobars  (solid  lines)  of  Fig¬ 
ure  10.  As  shown  by  the  streamlines,  the  plasma 
then  turns  back  so  that  it  is  flowing  clmcst  par¬ 
allel  to  the  anode.  This  turn  is  accompanied  by 
a  drop  in  Mach  number,  to  subsonic  speeds  along 
some  streamlines.  The  initial  acceleration  and  sub¬ 
sequent  drop  in  Mach  number  are  shown  in  Figure 
9,  contours  of  constant  Mach  number.  The  label  A 
is  at  a  local  maximum  of  1.57.  As  shown  by  the 
labelled  Mach  1  contour,  along  a  number  of  axial 
lines  the  plasma  is  dropping  from  supersonic  to  sub¬ 
sonic  speeds.  Along  streamlines  near  the  anode  the 
plasma  first  sees  decreasing  pressure  as  it  acceler¬ 
ates,  and  then  increasing  pressure  as  the  streamlines 
straighten  out.  The  current  lines  (solid  lines)  and 
constant  potential  contours  (dotted  lines)  shown  in 
Figure  II,  also  show  a  sharp  transition  in  this  region. 
Both  0.1  mm  &om  the  anode  and  along  a  ridge  of 
high  Hall  parameter  between  s  =  0.0015  and  s  = 
0.0075,  the  lines  run  almost  parallel  to  each  other, 
parallel  to  the  anode  very  near  the  anode,  and  par¬ 
allel  to  the  backplate  along  the  ridge.  Further  along 
the  channel,  as  the  pressure  increases  and  the  Hall 
parameter  decreases,  the  angle  between  the  lines  be¬ 
comes  larger. 

When  the  plasma  parameters  at  the  last  inte¬ 
rior  point  are  input  to  the  near  anode  boundary 
equations,  they  predict  significant  voltage  drops,  as 
shown  in  Figure  12,  transverse  cuts  of  the  potential 
drop. 

How  do  the  anode  voltage  drops  behave  at  dif¬ 
ferent  currents?  Heimerdinger[7}  directly  measured 
the  anode  voltage  drop  over  a  range  of  currents  in 
the  FFC  and  at  60kA  in  the  CAC.  The  measured 
voltage  drop  is  the  difference  between  the  potential 
measured  2  mm  from  the  anode  at  an  axial  loca¬ 
tion  of  0.043  m  &om  the  inlet,  and  the  potential  at 
the  anode.  His  data  are  shown  in  Figure  13.  The 
data  show  that  at  low  applied  currents,  anode  volt¬ 
age  drops  are  negligible  or  non-existent.  As  the  ap¬ 
plied  current  is  increased,  the  voltage  drops  appear 
and  increase  with  increasing  applied  current.  The 
anode  drops  seem  to  level  out  as  the  total  current 
is  increased  past  50-55  kA,  although  the  large  error 
bars  make  this  difficult  to  ascertain.  The  one  data 
point  available  for  the  CAC  seems  to  indicate  that 
the  anode  voltage  drop  in  this  thruster  is  somewhat 
smaller  than  in  the  FFC.  In  general,  the  anode  drops 
seem  to  account  for  50  -  75  %  of  the  total  terminal 
voltage  in  those  cases  for  which  the  anode  seems  to 
be  starved,  i.e.  for  currents  above  25  kA  or  so. 


The  numerically  predicted  voltage  drops  behave 
in  a  similar  manner.  Simulations  were  performed 
for  the  CAC  at  currents  of  23.4  kA,  27.3  kA,  31.2 
kA,  35.9  kA,  and  39.0  kA.  The  computed  numerical 
voltage  drops  are  also  plotted  in  Figure  13.  At  23.4 
kA,  the  voltage  drops  ate  quite  small.  As  current  is 
increased,  the  numerically  predicted  anode  voltage 
drops  grow,  like  the  experimental  data.  The  total 
potential  drops  from  the  experimental  data  and  the 
numerical  simulation  are  shown  in  Figure  14.  Agree¬ 
ment  gets  better  as  the  current  is  increased,  although 
whether  this  trend  will  continue  with  current  beyond 
39  kA  is  unclear.  The  reason  for  the  large  discrep¬ 
ancy  at  low  currents  is  also  unclear.  Possibly  this 
is  due  to  artificial  ignition  of  the  plasma  in  the  low 
current  numerical  cases,  due  to  the  lower  limit  on 
ionisation  fraction.  It  might  also  be  due  to  excessive 
damping  in  the  numerical  simulation  at  low  Mach 
numbers. 


The  analytical  results  demonstrate  that  starvation 
of  the  neat  anode  plasma  could  result  in  significant 
voltage  drops  in  the  quasi  neutral  plasma.  The  com¬ 
putational  results  show  that  this  mechanism  causes 
voltage  drops  with  behaviour  and  magnitude 
to  those  in  an  experimental  device.  As  these  voltage 
drops  are  a  major  cause  of  thruster  inefficiency,  un¬ 
derstanding  them  should  lead  to  improved  thruster 
designs. 
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bers  on  the  plot  show  values  at  the  local  minima 
(-582,-177,-61,-346  m/sec)  and  the  maximum  value 
(754  m/sec).  It  is  also  illustrated  by  the  stream¬ 
lines,  shown  as  the  dotted  lines  in  Figure  10.  This 
turn  leads  to  depletion  of  the  plasma  near  the  an¬ 
ode,  as  shown  by  the  isobars  (solid  lines)  of  Fig¬ 
ure  10.  As  shown  by  the  streamlines,  the  plasma 
then  turns  back  so  that  it  is  flowing  almost  par¬ 
allel  to  the  anode.  This  turn  Ls  accompanied  by 
a  drop  in  Mach  number,  to  subsonic  speeds  along 
some  streamlines.  The  initial  acceleration  and  sub¬ 
sequent  drop  in  Mach  number  are  shown  in  Figure 
9,  contours  of  constant  Mach  number.  The  label  A 
is  at  a  local  maximum  of  1.57.  As  shown  by  the 
labelled  Mach  1  contour,  along  a  number  of  axial 
lines  the  plasma  is  dropping  &om  supersonic  to  sub¬ 
sonic  speeds.  Along  streamlines  near  the  anode  the 
plasma  first  sees  decreasing  pressure  as  it  acceler¬ 
ates,  and  then  increasing  pressure  as  the  streamlines 
straighten  out.  The  current  lines  (solid  lines)  and 
constant  potential  contours  (dotted  lines)  shown  in 
Figure  11,  also  show  a  sharp  transition  in  this  region. 
Both  0.1  mm  from  the  anode  and  along  a  ridge  of 
high  Hall  parameter  between  z  =  0.0015  and  z  = 
0.0075,  the  lines  run  almost  parallel  to  each  other, 
parallel  to  the  anode  very  near  the  anode,  and  par¬ 
allel  to  the  backplate  along  the  ridge.  Further  along 
the  channel,  as  the  pressure  increases  and  the  Hail 
parameter  decreases,  the  angle  between  the  lines  be¬ 
comes  larger. 

When  the  plasma  parameters  at  the  last  inte¬ 
rior  point  are  input  to  the  near  anode  boundary 
equations,  they  predict  significant  voltage  drops,  as 
shown  in  Figure  12,  transverse  cuts  of  the  potential 
drop. 

How  do  the  anode  voltage  drops  behave  at  dif¬ 
ferent  currents?  Heimerdinger[7]  directly  measured 
the  anode  voltage  drop  over  a  range  of  currents  in 
the  FFC  and  at  60kA  in  the  CAC.  The  measured 
voltage  drop  is  the  difference  between  the  potential 
meitsured  2  mm  from  the  anode  at  an  axial  loca¬ 
tion  of  0.043  m  from  the  inlet,  and  the  potential  at 
the  anode.  His  data  are  shown  in  Figure  13.  The 
data  show  that  at  low  applied  currents,  anode  volt¬ 
age  drops  are  negligible  or  non-existent.  As  the  ap¬ 
plied  current  is  increased,  the  voltage  drops  appear 
and  increase  with  increasing  applied  current.  The 
anode  drops  seem  to  level  out  as  the  total  current 
is  increased  past  50-55  kA,  although  the  large  error 
bars  make  this  difficult  to  ascertain.  The  one  data 
point  available  for  the  CAC  seems  to  indicate  that 
the  anode  voltage  drop  in  this  thruster  is  somewhat 
smaller  than  in  the  FFC.  In  general,  the  anode  drops 
seem  to  account  for  50  -  75  %  of  the  total  terminal 
voltage  in  those  cases  for  which  the  anode  seems  to 
be  starved,  i.e.  for  currents  above  25  kA  or  so. 


The  numerically  predicted  voltage  drops  behave 
in  a  similar  manner.  Simulations  were  performed 
for  the  CAC  at  currents  of  23.4  kA,  27.3  kA,  31.2 
kA,  35.9  kA,  and  39.0  kA.  The  computed  numerical 
voltage  drops  are  also  plotted  in  Figure  13.  At  23.4 
kA,  the  voltage  drops  are  quite  small.  As  current  is 
increased,  the  numerically  predicted  anode  voltage 
drops  grow,  like  the  experimental  data.  The  total 
potential  drops  from  the  experimental  data  and  the 
numerical  simulation  are  shown  in  Figure  14.  Agree¬ 
ment  gets  better  as  the  current  is  increased,  although 
whether  this  trend  will  continue  with  current  beyond 
39  kA  is  unclear.  The  reason  for  the  large  discrep¬ 
ancy  at  low  currents  is  also  unclear.  Possibly  this 
is  due  to  artificial  ignition  of  the  plasma  in  the  low 
current  numerical  cases,  due  to  the  lower  limit  on 
ionization  fraction.  It  might  also  be  due  to  excessive 
damping  in  the  numerical  simulation  at  low  Mach 
numbers. 

6  Conclusions 

The  analytical  results  demonstrate  that  starvation 
of  the  near  anode  plasma  could  result  in  significant 
voltage  drops  in  the  quasi  neutral  plasma.  The  com¬ 
putational  results  show  that  this  mechanism  causes 
voltage  drops  with  behaviour  and  magnitude  similar 
to  those  in  an  experimental  device.  As  these  voltage 
drops  are  a  major  cause  of  thruster  inefficiency,  un- 
derstzmding  them  should  lead  to  improved  thruster 
designs. 
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Figure  9:  Mach  Number  Contours 
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Figure  12:  Radial  Cuts  of  Potential  Drop  in  31.2  kA 
CAC  (Baseline) 


Figure  10:  Stream  Lines  and  Constant  Pressure 
Contours 
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Figure  13:  Experimental  and  Numerical  Anode  Volt¬ 
age  Drops 
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Figure  11:  Current  Lines  and  Constant  Potential 
Contours 


Figure  14;  Experimental  and  Numerical  Total  Volt¬ 
age  Drops 
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Abstract 

For  purposes  of  systems  tradeoffs,  as 
well  as  for  preliminary  design  purposes,  tbere  is 
a  need  for  relatively  simple  methods  to  estimate 
aicjet  performance  wither  the  cost  and  delay  of 
a  full-fledged  numerical  study.  We  present  here 
one  such  method,  and  its  verification  against  a 
variety  of  test  data. 


NcmcnclatuM 

a  Ra/R 

Average  slope  of  ^  vs.  a 
£  Axial  electric  field 

F  Thrust 

C  Average  value  of  h  /  <r 

h  Static  enthalpy 

hg  Static  enthalpy  at  arc  edge 

H  Total  enthalpy 

H  Total  enthalpy  flux 

Hq  Enthalpy  flux  in  arc 

I  Arc  current 

J  Momentum  flux 

K  Thermal  conductivity 

m  Total  mass  flow  parameter  (Eq.  41) 
Arc  mass  flow 

N  Frozen  flow  patamer  (Eq.  4 1 ) 
p,  po  Pressure,  stagnation  pressure 

q  pu^ 

Q  Effective  value  of  y !  y-\ 

r  Radial  coordinate 

R  Radius  of  flow  passage 

Ra  Radius  of  arc 

T  Temperature 

u  Axial  velocity 

V  Arc  voltage 


^Professor,  member  AIAA 
‘Graduate  Student 


X  Axial  coordinate 

z  p/q 

y  Specific  heat  ratio 

e  h^ 

p  Heat  flux  potential 

p  Density 

a  Bectiical  conductivity 

1  Introducttoo 

The  model  is  based  on  a  two-stream 
approximation,  plus  a  set  of  simplified  reiadons 
for  gas  physical  propcrtes-eodialpy. 
conductivity,  etc.  Tlie  outer  stream  is  treated  as 
ideal  adiab^  flow,  with  Ohmic  heating  and 
heat  conduction  limited  to  the  inner  (”acc*) 
stream.  Growth  of  this  inner  stream  follows 
from  beat  diffusion  across  external  streamlines. 
By  posnilaiing  a  reasonable  ret  of  radial  shape 
functions  for  the  conductivity  mid  the  dynamic 

bead  (pu^)  the  conservation  laws  beocme  a  set 

of  (wdinary  differential  equadons,  which  must  be 
integrated  subject  to  prescribed  initial  obnditioos 
and  to  smooth  passage  through  the  sonic  sectian. 
Several  aspects  of  the  actnal  physia  aee  outtide 
the  scope  of  the  model,  and  most  be  externally 
supplied.  These  include  the  arc  attachment  point 
and  the  anode  voltage  drop,  plus  a  few  less 
significant  effects,  such  as  tte  chemical  fieeaing 
point  or  the  initial  cenieriine  eatha4ty.  'The 
results  presented  include  a  study  of  die 
sensitivity  of  predicted  perfbnumoe  to  these 
inputs. 

The  model  was  tested  against  an  extensive  set  of 
data  from  Stuttgart  University  on  a  2(Xcw 
thruster,  over  a  range  of  currents,  flow  rates, 
propellants  and  cooling  methods,  and  also 
against  data  from  two  higher  power  and  one 
smaller  power  thruster.  Results  indude  chamber 
pressure,  thrust,  voltage,  arc  radius  and  anode 
heat  loss,  plus  derived  quantities  such  as  specific 


impulse  aad  effideiicy.  The  agreement  witb  daa 
is  reasonable  good. 


I  ~E^\y2Krdr 


(7) 


2  Formulatiqn  and  Assumptions 

The  geometry  considered  is  shown  in 
Fig.  1.  The  region  outside  the  arc  (R>t>R^  is 
modeled  as  an  tsentropic  flow  region,  where  the 
gas  convected  from  upstream  has  not  yet  been 
affected  by  the  diffusive  spread  of  heat  coming 
from  the  arc  itself  (r<Ra). 


where  £  is  the  axial  electric  field,  and  O’  is  the 
electricai  conductivity  of  the  gas.  The  current  is 
kept  invariant  with  x  up  to  the  attachment  cross- 
section,  where  it  goes  abruptly  to  zero. 

The  wall  beat  flux  term  in  Eq.  (2)  will 
be  neglected,  except  in  the  anode  attachment 
region  as  will  be  discussed  below.  Accordingly, 
Eq.  (9)  gives  for  the  electric  potential  from  the 
cathode  tip 


This  is  supported  by  photometric  data  ^  ^  ^  which 
show  a  sharply  delineated  arc  boundary 
surrounding  by  non-emitting  gas,  and  also  by 
results  of  detailed  2-D  simulations  according 
to  which  the  buffer  gas  is  actually  colder  than  the 
walls,  with  a  thin  w^l  boundary  layer  of  the 
order  of  0.3  mm  total  thickness  (for  a  constrictor 
diameter  of  2.5  mm). 

The  arc  itself  is  a  region  of  intense 
Ohmic  heating  and  radial  heat  conduction. 
Radiant  heat  transfer  is  comparatively  negligible 
for  arc  diameters  less  than  about  1  cm. 
Convection  is,  of  course,  the  other  important 
effect.  Ignoring  friction,  the  governing  equations 
(conservation  of  mass,  energy  and  momenuim) 
can  be  integrated  across  the  flow  area  (o<r<R(x)) 
to  obtain 


where  the  overall  fluxes  m,  H  and /an 


n 

m  *  \qPu  Ixrdr 

(4) 

pH  I  h  +  I  litrdr 

(5) 

J  =  +  2ftrdr 

(6) 

The  total  current  1  in  Eq.  (2 

)  is 

calculated  from  Ohm's  law  as 

V-^Vc 


H-H{o) 

I 


(8 


where  AVc  is  the  cathode  potendal  drop,  and  H 
(0)  is  the  axial  enthalpy  flux  at  the  cathode  tip 
section.  In  the  anode  j-tiacfainent  region,  the 
quasi-one  dimensionality  is  locally  violated,  and 
a  large  radial  electric  field  does  occur  leading  to 
the  observed  anode  drop.  AVa.  Hie  energy  AVa. 
dissipated  in  this  region  is  mainly  absorbed  by 
the  wall,  and  tbe  axial  enibalpy  flux  H  must 
therefore  by  -AAT  •  I AV^. 


The  axial  evobirioo  of  the  vc  mass  flux 
ffig  which  is  deflned  as  in  Eq.  4.  but  with  the 
upper  limit  changed  to  Ra(x),  ia  governed  by  the 
diffusive  spread  of  the  region  of  heated  gas, 
which  covers  a  steadily  increasing  number  of 
stream  tubea.  This  bet  will  be  used  bere  to 
formulate  a  simpie  growth  model  Eveoibougb 
the  enthalpy  radial  profile  h(r)  is  in  reality 
continuous,  it  is,  as  we  noted,  very  sieqi  at  tbe 
edge  of  the  arc.  This  partly  due  to  tbe 
counterflow  heat  diffhaiOB,  but  abo  the  very 
lapid  variuriaa  of  electrical  condnetivity  whh 
temperature  in  the  5000-7000  Krqkm,  where 
the  gas  begins  to  iraiize  tbennaily.  Since  tbe 

beating  rate  is  aiT)E  ooceatenqieratnre  r«~- 
6000  K  is  reached,  dissipatioo  increases  rqpidly 
with  further  T  increases,  and  the  T-profile  must 
steepen.  Consider  an  idealized  situation  as 
shown  in  Fig.  4,  where  he  (x)  is  the  arc  centerline 
enthalpy,  hout  U)  is  the  adiabatic  outer  entba^ty, 
and  he,  cotiesponding  to  Tg  36000  K,  is  the  arc- 
edge  enthalpy.  If  the  arc  is  "ingesiing"  gas  at  the 


dih/j  /  \  diha 

rate  — tbe  amount  (rf* heat  ~^out) - 

cbi  dx 


must  be  supplied  to  this  region  pa  unit  length 
(more  precisely,  toigl  enthalpies  should  be  used. 


including  kinetic  energies).  This  heat  is  supplied 
by  radial  cooducdoo.  If 


is  (he  radial  heat  flux  at  the  arc  edge,  we  obtain 

(^e  -  ^out ) =  2;r  Ra  Qc  (9) 

Turning  now  to  the  momentum 
equation,  the  overall  conservation  law  (Eq.  3) 
needs  to  be  supplemented  by  some  other 
equation  which  expresses  the  character  of  the  arc 
axial  momentum  balance.  Since  the  dynamic 
2 

bead  pu  which  develops  gradually  in  the  arc- 
heated  duct  is  generated  by  the  gradual  drop  in 
pressure  and  p  is  radially  uniform,  it  can  be 

2 

reasonably  expected  that  pu  itself  will  also  be 
nearly  independent  of  r.  This  is  confirmed  by 
inspection  of  full  2-D  solutions,  such  as  those 
shown  in  Refs.  2.4. 


In  the  interest  of  algebraic  simplicity, 
we  would  like  to  interrelate  the  gas  enthalpy, 
thermal  and  electrical  conductivity,  pressure  and 
density  in  as  simple  a  fashion  as  possible,  while 
preserving  the  b^ic  physics.  The  following 
approximations  have  bwn  found  reasonable  and 
useful: 

(a)  h  =  Q-  (10) 

P 

y 

Of  course,  for  an  ideal  gas,  Q  = - .  y  being 

y-1 

the  specific  heat  rado.  F(»  a  diatomic  gas,  then, 
Qs3.S,  and  for  a  monatomic  gas  Qs2.5.  Fora 
non-ideal  gas  undergoing  dissociation  and 
ionization,  /  will  approach  unity,  and  Q  can  be 
rather  large  over  limited  T  intervals.  From  data 
collected,  for  example  by  Lee  (Ref.  S)  Q  is  on 
average  7-8  for  H2  and  N2  a  broad  interval, 
to  over  30.0(X)K.  The  parameter  Q  (and  others) 
is  given  in  Table  1  for  N2  and  H2  .  as  well  as  for 
N2H4,  for  which  simple  molar  averaging  has 
been  used. 


(b)  The  tbennal  conductivity  K  has  very 
sharp  peaks  (due  to  molecular  diffusion  followed 
by  recombination)  in  the  dissociation  and 
ionizadon  ranges.  Fortunately,  the  quantity  of 
interest  is  the  so-called  heat-flux  potential. 

which  varies  more  smoothly.  In  fact,  p  is 
almost  linearly  related  to  the  electrical 
conducdvtty  a.  as  shown  in  Fig.  3  for  Hydrogen 
at  I  atm.  We  adopt  the  form 

^  5  (12) 

with  parameters  and  B  as  shown  in  Table  I . 

(c)  It  is  also  useful  to  relate  the  enthalpy  b 
to  the  electrical  conducdvity  p.  Fig.  (4)  shows 
this  relationship  for  hydrogen  at  I  atm.  To  a 
rough  approximadoo,  we  will  use 

hsGff  (13) 

with  G  also  ^lecified  in  Table  1 

These  coneladons  are  reladvely 
insensidve  to  pressure;  some  correction  is 
included  in  the  model  to  account  for  the  Spitzer 
logarithm  effects 

The  radial  variadons  of  velocity,  mass 
flux,  enthalpy  flux  and  kinetic  energy  can  now 
be  expressed  in  terms  of  only  those  of  the 
enthalpy.  We  define 


9(i)*pu^ 
and  obtain 


(independrat  of  r)  (14) 


(15) 


puh  =  yjOpq  Va 


2  2  2  \Qp 

and  hence 


puhj  s  pul  ^  ~ 


(17) 

(18) 

(19) 


In  these  expressions,  the  factors 
involving  Q,  p  and  q  are  independent  of  r.  These 
expressions  indicate  a  large  velocity,  but  a  small 
mass  flux  inside  the  arc.  where  h  is  high.  The 
enthalpy  flux  is  large,  but  only  in  proportion  to 

Va  rather  than  to  h ,  as  one  might  have 
expected. 

It  is  clear  then,  that  the  speciftcation  of 
the  radial  variation  of,  say.  electrical 
conductivity,  would  suffice  to  specify  those  of  ail 
the  quantities  needed  to  evaluate  the  integrals 
rh,  m^,  H  (since  (T  is  related  to  h  by  Eq.  (13). 
A  two-parameter  profile  is  adopted,  of  the  form. 


(20) 


and  a  convenient  form  for  f  is  obtained  by 
consideration  of  the  simpler  problem  of  an  arc 
which  fills  a  cylindrical  cooled  wall  of  radius  Ra: 


not  a  good  approximatioo.  although  it  does  show 
the  correct  trends. 

Because  of  Eq.  (13).  this  also  gives  the  enthalpy 
profile; 

A(r.j)a  Aj.(x)  /of  2-405--^  I  (22) 


Substituting  Eq.  (21)  into  the  total  current 
expression  (Eq.  1)  gives 


I^IkE  Oc  Jo  2.405—  1  rdr 


Rn 


=  l.3S6EacRa 


(23) 


This  expresses  the  axial  field  E  and  the 
dissipation  per  unit  length,  El,  in  terms  of  I, 
(T^and/fg.  The  heat  flux  at  the  arc  edge 
needed  in  Eq.  (9)  can  also  be  evaluated  now; 


or  licR^qf.  =  7.843 


(24) 


Finally,  we  can  use  Eq.  (22)  to  evahiaK 
the  mass  and  energy  flux  integi^  the  arc 
mass  flow. 


ffig  =  /o®pH  2«r  dr  = 


2>r 


=  7. 634 Vo 


(25) 


(T  =  (T^Jq 


2.405 — 


(21) 


where  Jo(y)  is  the  zero'th  order  Bessel  function, 
whose  first  zero  is  at  ys2.405.  In  the  purely 

B 

cylindrical  case,  the  relationship  E=2.405  —  is 

Rq 

obtained  from  the  radial  heat  balance.  Since  we 
now  have  substantial  convective  cooling,  this  is 


The  mass  flow  rate  outside  the  arc  is 
easier  to  evaluate,  because  hahout^cotur.  in  that 
region.  We  obtain 

ffi  -  flig  =  ipu)^ut  -  R^)  » 

_  R^-R^ 

ffV^VW-T.'— (26) 


m  »  const. 


(29) 


wbere  we  have  used  Qout  of  Q  value, 
which  is  appropriate  only  to  the  high  tempeiature 
region.  For  A/j  or  we  can  use  Qout*3.5 
( y-1.4). 

For  the  total  enthalpy  flux  in  the  arc. 


X  tix] 


the  new  integral  is  calculated  to  be  1 .762.  giving 


Again,  the  enthalpy  flux  outside  the  arc 
is  simpler.  Using  Eq.  (19)  we  obtain 


We  compile  here  the  relevant  governing 
equations.  Non-dtmensional  arc  radius  and 

enthalpies  are  introduced  as  a  = 

R 

9^=^  ,  9^-  and,  using  Eqs. 

(1,2)  (with  zero  wall  heat  flux).  (3)  and  (9),  and 
expressing  the  dissipadon  and  the  arc-edge  beat 
flux  fixMn  Eqs.  (23)  and  (24),  we  obtain 


dm  ,  -  ,  0 

—  =  7.843 - 

dx  Cl-  001^1 

(30) 

dJ  dR 

—  =  p2xR  — 

(31) 

dx  dx 

dH  C  P-  1 

(32) 

- - - - - - 

dx  1.356  hg  df.  R'^ 

®Ottt  ~  ®oux,o  I 
\PoJ 

The  integral  quandties 

(m,  JassiH)  are  related  to  the  primidve 
variables  (p. q.  a.  9^) by 

J  Mtf  R^{p  +  q)  (36) 

Thus,  after  advancing  one  step  in  the 
axial  direcdon.  the  set  of  non-linear  equadoos 
(34-37)  must  be  solved  for  (p.  q.  a.  0g)wbkk 
then  can  be  used  to  evaluate  the  axial  derivatives 
and  repeat  the  cycle. 

Equadoos  (34-37)  can  be  reduced  to  a 
single  nonlinear  equadao  in  the  quantity  Si*q/p. 
which  for  an  ideal  gas  would  be  equivalent  to 

,M  being  the  Mach  number.  The  quantities 
/hg,  /f  and  7.  which  are  knowa  from  the 

integradon's  to  the  current  x  locatioo,  are 
parameten  in  this  equation  for  z.  This  equstion 
can  have  two,  one  or  no  roots.  Starting  from  the 
low-speed  rod,  wbere  is  small  two  tooa 


occur,  of  whidi  ooe  can  be  identified  as  the 
desired  subsonic  solution.  As  increases,  and 
if  is  too  higii  for  tbe  givoi  m.  tbe  end  of  tbe 
constrictor  is  reached  within  this  subsonic 
branch.  Continuation  into  the  widening  nozzle 
region  then  leads  to  a  Venturi-type  solution, 
which  fails  to  expand  tbe  flow  supenonically. 

On  tbe  other  hand,  a  Po  which  is  too  low  leads  to 
a  situation  where  no  solutions  exist  beyond  a 
certain  point  in  tbe  constrictor.  Tbe  correct  is 
such  as  to  place  tbe  sonic  point  at  the  constrictor 
exit,  such  that  tbe  sudden  area  divergence  allows 
Che  solution  to  transition  smoothly  into  the 
supersonic  branch.  The  solution  itself  is 
accomplished  at  each  x  value  in  a  few  Newton- 
Raphson  iterations. 

Once  zaq/p  is  known,  the  remaining  primitive 
variables  can  be  easily  found; 

J  s 

p  *  q  m  -  —  J  ■  i  ;  p  ■  — —  ;  q  m  zp  (38) 

x  K  1  +  r 


can  be  shown  that  this  effea  can  be 
characterized  by  tbe  parameter 


(41 


which  reduces  to  unity  when  Qi^aQ.  Tbe  factor 

I  +  ““  in  Eq.  (27)  (for  should  now  be 
2Q 

alp 

replaced  by  1/  N  — N,  while  a  bctor  1/N 
2Q 

should  multiply  tbe  right-hand  side  of  Eq.  (25) 
(for  fhj).  Similarly,  Eq.  (39)  should  have  a 

factor  of  l/N^,  and  Eq.  (40)  a  factor  •Jn  .  The 
iteratioo  for  z  is  now  slightly  more  involved, 
because  N  itself  depends  on  z.  but  this  is  easily 
accounted  for  by  using  Eqs.  (39)  and  (41)  at  each 
step  of  the  iteration. 


m  3.493 


-i2 


0.4266  ^ 

Vz  i/C/\  y  /  m 


(39) 


(40) 


out 


4  Nogglt  Eapanslon.  Ftoan  flaw  and 
Other  Effwts 


The  physics  f  the  anode  region  is 
beyond  the  s<^  of  this  model  Experimeniil 
evidence  suggests  that  the  arc  anachei  oev  the 
beginning  of  the  supersonic  noaaJe  if  the  walls 
ate  hot,  and  somewhat  (uiher  downstream  if 
theyaiecold.  The  voltage  drop  AVy|  in  this 
attachment  is  also  difScuU  10  predia  ftom  first 
princ^les.  Experimentt  on  the  Gennan  TTl 
thruster  (2]  which  is  waier-GOOled,  and  on  its 
radiation  cooled  counte^pmtt^l  show  that  AVa 
is  higher  by -'20V  in  the  case  wheiB  the  walls  ate 
cold. 


Once  tbe  flow  enters  tbe  divergem 
nozzle,  tbe  pressure  falls  rapidly,  and  chcmiral 
reaction  rates  fails  to  keep  up  with  the  evolving 
conditions  (frozen  flow).  From  the 
thermodynamic  viewpoint,  tbe  principal 
consequence  of  this  is  that  the  failure  to 
recombine  keeps  large  amounts  of  beat  of 
recombinadon  from  becoming  thermal  energy, 
which  could  then  be  converted  to  directed  kinetic 
energy  by  tbe  nozzle.  Formally  speaking,  the 
curve  of  enthalpy  h  vs.  p/  p  will  have  a  slope 

s  Y  I{y-\)  corresponding  to  a  non- 

reacting  mixture  of  mono-  and  diatomic 
molecules,  rather  than  the  much  steeper  Q  which 
is  used  for  tbe  hot.  reacting  gas  in  tbe  arc 
upstream  of  the  freezing  point 

Assuming  h  is  still  linear  witb  p/  p ,  but 
with  tbe  new  slope,  and  further,  that  fluid 
elements  move  ^ong  lines  of  constant  r/Ra(x).  it 


present  model  shows  that  the  fracdon 

condnoes  to  increax  in  the  noak, « it  should  in 
the  presence  of  heat  diffiisioo.  This.pha(ihB 
expansion  of  R(x)  itself,  makes  the  'arc”  radius 
increase  rapidly,  reducing  its  cooling  rate,  and 
hmice  its  axial  electric  field.  Beyoodefew 
constrictor  diameters  iiuo  the  nozzle,  (he  sic 
voltage  increases  very  little,  aad  arbitrarily 
moving  tbe  attachment  point  iq>  or  downstream 
has  only  a  minor  effect  on  total  volo^.  Fbribe 
numerical  appikarions  reported  here,  we  heve 
assumed  attachment  at  the  constrictor  exit  for 
bot-wall  cases,  and  about  2  constrictor  diagieiers 
into  the  nozzle  for  cold  walls. 

Wall  fricdon  is  a  relatively  small  effect 
in  arcjets.  but  it  does  have  some  effea  on  tbe  na 
peribnnanoe.  maiiily  by  reducing  the  nozzle 
efficiency.  We  have  used  a  relsrively  crude 
model  to  account  for  it.  This  consists  of  adding 


to  tlM  integral  momentum  equation  a  frictional 
term  derived  from  tbe  numerical  simulations  of 
Ref.  3.  For  tbe  constrictor,  ibis  amounts  to  a 
friction  factor  ^iO.OlS.  In  tbe  nozzle,  c/ is 
allowed  to  decrease  according  to  a  flat-plate 
formula . 


S  Arciet  Performance.  Method  of 

Cakutatifln 

The  vacuum  thrust  of  the  arqet  is  given 

by 

F  =  \oexa  piP"  Inrdr  +  {pnR^ (42) 


and,  in  fact  tbe  impulse  I(x)  gives  the  thrust  that 
would  result  by  terminating  the  nozzle  at  the 
current  x  station.  SpeciTtc  impulse  and  efficiency 
then  follow  from  standard  definitions. 
Computations  are  performed  according  to  tbe 
following  sequence: 


( 1)  Specify  geometric  parameters,  nozzle 
contour,  etc. 

(2)  Specify  initial  outer  temperaure  Tgutfo) 
(acmally  at  the  cathode  tqi). 

(3)  Prescribe  the  anode  attactoent  point 
Xatt  and  tbe  anode  voltage  drop  AV^. 

(4)  Specify  mass  flow  m  and  current  I. 

(5)  Guess  an  initial  pressure  Po- 

(6)  Integrate  forward  to  tbe  constktor  exit 
and  check  for  smooth  sonic  passage.  If 
not  achieved,  modify  po  and  repeat. 

(7)  When  smooth  sonic  passage  achieved, 
continue  integration  to  attachment  pome 
At  xaxait,  add  dV/^.  to  voltage,  subtract 
LIVa  from  H. 

(8)  Contiaae  to  nozzle  exit.  Calculate 
thrual.  specific  impulse 


1  F 

lj„  - - ,  efficiency  ij  = 

gm 


2m  IV 


Some  additional  details  must  be 
discussed  about  tbe  initial  values  used  to  start  tbe 
computation.  As  noted,  tbe  chamber  pressure  Po 
is  found  by  the  sonic  passage  condition.  Less 
clearly  defined  are  the  arc  parameters  (radius  and 
temperature)  right  at  tbe  c^ode  tip  where  the 
integration  starts.  It  is  found  that  a  fairly  wide 


latitude  exists  in  tbe  choice  of  a(o)«Rg(oVR  and 
of  6f.{o)  =  H^{o)  /  hf .  The  effect  of  a  change  in 

either  of  these  parameters  is  only  felt  in  tbe  near 
vicinity  of  tbe  cathode,  and  tbe  arc  behavior  in 
tbe  constrictor  is  very  nearly  independent  of 
these  changes.  Tbe  procedure  adopted  is 
therefore  to  select  a(o)  and  &c(o)  such  as  to 
generate  smooth  profiles  of  a(x),  (j:)  near  tbe 

cathode.  Typically  this  leads  to  a(o)~0.2-0.3  and 
0c(jt)-15-3O  (basedonheSh(6000K)forthe 
gas  considered). 


6  Results  and  Comparison  tn  Data 

Tbe  medium-power  (6-30  kW)  arqeu 
of  References  (1,2,6)  have  been  very  well 
documented  experimentally,  and  constitute  an 
excellent  test  bed  for  theoretical  models.  They 
cover  a  wide  range  of  flow  rates  (0.05-0.3  g/sec), 
currents  (50-200  Amp)  propellaaa  (Hz,  Nz  * 

Hz,  Argon)  and  cooling  arrangements  (water 
cooling,  radiation  cooling).  Measured  quaotitiet 
include  voltage,  stagnation  picswie,  beat  kiu 
and  its  distribution,  thrust  and,  in  some  cases,  arc 
radius  and  constrictor  ptessme. 

The  basic  TTl  thruster  is  water-cooled 
with  segmented  copper  nozzle  and  constrictor. 
Tbe  constrictor  has  a  length  of  Smm  and  a 
diameter  of  2.5mffl.  Tbe  area  ratio  of  the 
supersonic  nozzle  is  lOfkl  with  a  17.5*  mean 
nozzle  angle.  The  cathode  has  a  conical  tungsten 
tip  with  a  30*  cone  angle,  nuting  to  a  similar 
conical  chamber  wall  from  which  it  is  separated 
by  a  2  mm  axial  gap. 

Caicaiations  and  data  for  the  extreme 
values  of  mass  flow  rate  and  current  reported  in 
Ref.  [2]  ate  compiled  in  Table  2.  As  noted  in  tbe 
table,  attachment  was  throughout  assumed  at  xw9 
mm  (3.2  mm  past  tbe  constrictor  end).  The 
anode  drop  dV^s28V  was  kept  the  same  for  all 
cases,  and  was  selected  so  as  to  provide  tbe  best 
match  to  measured  voltages.  As  tbe  Qt 
column  shows,  this  AV/{ ,  times  tbe  total  ctarrent, 
also  matches  reasonably  well  against  the 
experimentally  measured  total  heat  lott.  The 
selected  values  of  bc(o)  and  a(o)  are  arbitrary 
(see  discussioo  in  Sec.  5)  but  insensitive. 

Tbe  initial  outer  temperatures 
Tout.o*7(X)  K  at  0. 1  g.sec,  500K  at  0.3  g/aec)  are 
reasonable  values  for  tbe  gas  which  has  been  in 
contact  with  (on  the  outside)  the  water-cooled 
copper  wall,  and  (on  the  inside)  the  hot  tungsten 


cathode.  It  was  found  that  the  chamber  pressure 

Po  varied  with  Tooto.  roughly  as  This 

is  due  to  the  fact  that  most  of  the  mass  flow  at 
Che  cbokiag  plane  is  in  fact  in  the  outer  flow 
region,  so  that  by  the  choking  relation 

m  *  PqA  I  •  Po  *  'j^OUt.O  * 

fixed  m.  Because  of  this  the  thrust  F  also 
increased  with  Tout.o~  The  good  match  achieved 
for  both  Po  and  F  against  experiment  validates 
the  selecticM)  of  Tout.o-  The  higher  Tout.o 
selected  for  the  lower  flow  rates  attempts  to 
reflect  the  higher  (power/mass  flow)  for  such 
cases,  and  does  appear  to  yield  proper  results. 

Inspection  of  Table  2  shows  that  all  the 
performance  parameters  are  predicted  to  within 
10%  (most  of  them  significantly  better  than  this). 
Of  course,  the  voltage  match  is  largely  enforced 
by  Che  choice,  and  that  of  the  thrust  is 
helped  by  the  rou/,0.choice.  However,  the  latter 
is  only  moderately  affected  by  variations  of 
Tout.o-  within  the  physically  reasonable  range. 
For  example,  for 

m  -  0.3  g  /  sec  and  I  «  ISO  A.  using 
Tout.o*^^^  instead  of  the  500K  of  Table  2 
would  change  the  thrust  from  1.64  to  1.80.  a  10% 
variation. 

Trends  with  mass  flow  and  current  are 
well  reproduced.  In  particular,  the  slight 
decrease  of  voltage  with  current  is  accurately 
predicted.  The  increase  in  voltage  with  mass 
flow  is  also  well  reproduced. 

Ref.  (I)  reported  experiments  in  which  a 
2mm  diameter  observation  window,  centered  2.6 
mm  downstream  of  the  cathode  tip  was  used  to 
optically  observe  the  arc  and  also  to  measure 
static  pressures.  A  comparison  of  measured  and 
calculated  arc  widths  is  shown  in  Fig.  5.  The 
excellent  agreement  further  verifies  the  model's 
validity. 

For  the  radiation-cooled  case,  repotted 
in  Ref.  (1),  the  geometry  was  the  same,  but  a 
monolytbic  tungsten  anode  was  used  instead  of 
the  water-cooled  segments.  The  outside 
temperature  of  the  anode  block  was  measured 
and  varied  from  1200'C  to  18(X)‘C  depending  on 
cunent.  Internal  temperatures  were  modeled 
numerically  and  were  estimated  to  yield 
Tgas  =  800“  C  at  the  chamber  inlet  a  low  power 

(2kW)  case.  This  is  not  suflicient  information  to 
calculate  Tout.o  accurately,  but  a  reasonable 


estimate  can  be  made  of  Tout,o  *  I400K  for 
n-O.l  g/scc. 

The  anode  beat  loss  was  not  directly 
measured  in  these  ussts.  but  the  observed  voltage 
decrease  of  30-40V  with  respect  to  the  water- 
cooled  case  suggests  that  the  anode  drop  must 
have  become  very  small  indeed.  In  fact  even 
allowing  for  >o  and  for  an  earlier  arc 
attachment  (at  the  constrictor  exit),  we  cannot 
obtain  such  a  large  decrease  in  voltage,  as  shown 
in  Table  3.  As  a  consequence,  even  though  the 
thrust  prediction  is  fairly  accurate,  the  efficiency 
is  underestimated  by  the  model.  Tbechamb^ 
pressure  was  not  measured,  and  it  is  conceivable 
flku  Touto  should  be  raised  somewhat  in  order  to 
effea  a  better  thrust  predicdon  by  increasing  po. 
but  this  would  exacerbate  the  voltage 
overpredictioa  problem.  It  is  possible  that  the 
anode  attachment  may  occur  even  before  the  end 
of  the  constrktor,  or  that  at  least  a  pan  of  the 
current  is  teaching  the  anode  there;  however,  the 
data  were  insuffldent  to  verify  this. 

Some  calculations  were  also  done  for  a 
siinilar  radiatioo -cooled  30  kW  bydrogeo 
thrusta  Grom  NASA-Lewis  RC  (Ref.  8).  This 
had  a  coosiiictor  radius  of  0.89  nan,  coQstrictior 
length  of  3.S6oun  and  a  20*  nozzle  with  270:1 
expansioa.  Rmulu  are  shown  in  Table  4.  Once 
again.  is  required,  but  here,  both  ihrnst 
and  voltage  are  weO  predicted  once  this  choice  is 
made. 

The  results  for  the  higher-powered 
HIPARC  tfartisaerPJ»  which  had  noonsnictor  of 
3  mm  ladias  ad  6  mm  length;  wilfe  a  1 14: 1, 20* 
nozzle,  were  similarfy  accuraie  aiO.2  g/sec,  3(X) 
A  :Vpredictedw94.6V  (includhig  a6V  dV^)  vs. 
Vmcasuieda96V,  aDdF|jnd.*l-80M  vs, 

Pn,eaf  .>I  .96N.  Oi  the  Other  hand.  St  0.3  g/sec, 
800A..  while  the  same  AVa  leads  m 
Vpfed.sl02V  vs.  Vmffw  »103V.  the  model  gives 
Ppied.«2.64N,  while  the  data  of  Ref.  (^1  intBcaiei 
Fmfu  «3.4(M4.  Inclusion  of  self-magnetic 
(MFD)  effects  acoountt  for  about  0.2N  only,  so  a 
substantial  disctepaacy  remaiBS.  It  should  be 
noted  tbat  the  tbrust  coefficient  implied  by  the 
data  is  (CF)aKa(.*2.46.  much  larger  than  the 
values  seen  in  other  thrusters  and  in  the  model, 
and  perhaps  unrealistic. 

Even  though  the  molar  averaging 
procedure  used  u  generate  the  coefficients  for 
N2H4  in  Table  1  may  be  open  to  criliciam,  the 
model  is  able  to  make  accurate  predkiions  for 
hydrazine-fueled  tbrusten  as  well  This  is  shown 


in  Table  5,  which  pertains  again  to  the  TTl 
water-cooled  thro^  of  Ref.  [1].  As  in  the 
voltage  data,  in  this  case  a  constant  20V  for  all 
cases.  Chamber  pressures  and  tlmists  are  then 
well  predicted  for  a  range  of  Qows  and  currents. 

7  Condusiona 

A  fairly  simple  model  which  retains 
only  the  essential  facts  about  the  physics  of  the 
arc  and  the  flow  process  has  been  shown  to  be 
able  to  predict  most  of  the  important 
performance  parameters  of  arcjets  with  an 
accuracy  comparable  with  water-cooled  H2  case, 
a  substantial  anode  drop  is  seen  to  be  required  to 
match  the  refined  numerical  simulations.  A 
notable  exception  is  the  model's  inability  to 
calculate  anode  drops,  for  which  the  near-anode 
region  needs  to  be  modeled  to  the  extent  of 
including  temperature  non-equilibrium  and 
charge  carrier  diffusion  (Ref.  3).  A  second 
shortcoming  is  the  need  to  provide  separate 
estimates  of  the  temperature  of  the  gas  outside 
the  arc;  this  should  be  amenable  to  relatively 
simple  thermal  analysis,  however.  The  model 
could  prove  useful  for  preliminary  design  and  for 
broad  systems  studies,  where  the  high  costs  of 
detailed  simulations  may  be  still  prohibitive. 
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Table  1:  Fit  approximations  for  hydrogen,  hydrazine  and  nitrogen 


G  (J  m/s  kg) 


2.47  X  103 


1.73  X  Vfi 


120x10^ 


Table  2:  Contputanonal  and  experimentti  data  comparison  for  watg-cooled 
20  kW  levd  thruster  with  hydrogen  propellictt 


Computaiioael 

(Experimentti) 

Case  lA 

(Duels 

Case  1C 

Case  ID 

m  (kg/sec  •  10-3)  0.1 


I  (amperes) 

Toafco(K) 
hco  (J/kg) 
Po(N/m2) 

V  (volts) 
F(N) 
l8o(«C) 

n 

QlM(kW) 

R./R@ 

consihcunrexit 


SO 

700 

SxlO’ 

9.02  xlO< 
(9.00  X  10*) 

165 

(165) 

0.58 

(0.60) 

592 

(612) 

0.204 

(0.218) 

1.40 

(1.30) 

0.598 

1.309 

(1.358) 


0.1 

150 

700 

8xlOP 

1.23  xlOS 
(IJOxlO*) 

155 

(140) 

0.81 

(0.83) 

876 

(847) 

0.141 

(0.164) 

4.20 

(4.00) 

0.758 

1.342 

(1.301) 


3X10P 

8xlOP 

2,03x10’ 

144x10’ 

(235x10’) 

a74xl0’) 

225 

199 

(233) 

(198) 

132 

1.63 

(1.17) 

(1.62) 

449 

554 

(398) 

(551) 

0358 

ai48 

(0.196) 

(P.147) 

1.40 

430 

(1.00) 

(330) 

0.462 

0.610 

1324 

1361 

(1.002) 

(1.186) 

For  all  cases.  AV  •  28  V,  Rx(0)  ■  02  nun,  xm  «  9  am,  Q  ■  7,  Qott  *  3  J,  md  Qfr  •  4.7S. 


Table  3:ConipKitatioaalandexperinKntalcUacotnpahso 
cooled  2(UW  level  thruster  with  hydrogen  propeUant 


Gmtputatioaal 

(Experimental) 

Caae2A 

CaaeZB 

m  (k^aec‘10’3) 

0.1 

0.1 

((aiQperes) 

50 

150 

ToutoCK) 

1400 

1400 

hCoC^W 

SxlOP 

SxlOP 

Po(NAn2) 

U7xl05 

1.31  X  lO^ 

V(vola) 

14S.2 

(120) 

WJ 

im 

F(N) 

0.73 

(0.73) 

0.96 

(0.99) 

743 

(763) 

979 

(1009) 

n 

0.367 

(0.467) 

0,261 

(0J03) 

Ra/R^ 

oonurictorexit 

0.538 

0.722 

Q* 

1.267 

1,299 

F0r«U cases,  AV.  •  0 V. RgCO) ■  OJ «». xan •  6 nun. Q - 7. Qa*! •  33. Q|r»4.7S. 


Table  4:  Computanonal  and  expenmenal  data  compazison  for  ladianon 
cooled  3(XJcW  level  thruster  with  hydrogen  propellant 


Congnitaiioiiai 

(Experimental) 

Caae3A 

Caae3B 

m  (kg/sec*  10*3) 

0.1228 

0.1228 

Kan^eres) 

107.x 

2603 

TooijodC) 

2000 

2000 

(Kod/kg) 

6x10^ 

9x10^ 

Po(Wn?) 

3.63x10$ 

537x10$ 

V(yohs) 

184 

(172) 

164 

(161) 

F(N) 

IJ238 

(1.238) 

1.70 

(1.761) 

I«(sec) 

1028 

(1028) 

1412 

(1461) 

7 

0331 

(Q331) 

0375 

(0398) 

Cf 

1.37 

135 

ForaUcaaes.  AV.  «OV.R«(0)sOJiniB,xatt«6inm.Q*7.QiMi«3J.aii 


Fig.  3  Calculated  heat  flux  potential 
vs.  electrical  conductivity 
(hydrogen,  p*l  atm.) 


